Figure Plotting Rev1

Notebook to make the panels for the paper figures after revision

General Setup

Setup chunk

Load libraries

knitr::opts_chunk$set(fig.width = 8)
knitr::opts_knit$set(root.dir = normalizePath(".."))
knitr::opts_knit$get("root.dir")

[1] "/nas/groups/treutlein/USERS/tomasgomes/projects/liver_regen"

Define plot theme(s)

library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(Seurat)

Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio
Registered S3 method overwritten by 'spatstat.geom':
  method     from
  print.boxx cli 
Attaching SeuratObject

library(ggplot2)
library(ggrepel)
library(ggridges)
library(pheatmap)
library(RColorBrewer)
library(MetBrewer)

Useful functions

th_gen = theme(axis.text = element_text(size = 7, colour = "black"),
               axis.title = element_text(size = 8, colour = "black"),
               axis.line = element_line(colour = "black"),
               plot.background = element_blank(),
               axis.ticks = element_blank(), 
               panel.background = element_blank(),
               panel.grid = element_blank(),
               legend.background = element_blank(), 
               legend.key.size = unit(0.35, "cm"),
               legend.key = element_blank(), legend.spacing = unit(0.7, "cm"),
               legend.text = element_text(size = 9),
               legend.title = element_text(size = 10),
               legend.box.margin = margin(0.1,0.1,0.1,0.1),
               legend.box.spacing = unit(0.5,"cm"))
pointsize = 0.3
theme_set(th_gen)

Colour palettes

breakStr = function(s, n = 20) {return(gsub(paste0('(.{1,',as.character(n),'})(\\s|$)'), '\\1\n', s))}

getTopTerms = function(godf, topt = 100, ncl = 5, nt = 2){
  topt = if(nrow(godf)<topt) nrow(godf) else topt
  if(nrow(godf)<ncl) return(godf)
  df = godf[1:topt,]
  genes = sapply(df$geneID, function(x) strsplit(x, "/"))
  resmat = matrix(0, length(genes), length(genes))
  for(i in 1:length(genes)){
    for(j in 1:length(genes)){
      resmat[i,j] = length(intersect(genes[[i]], genes[[j]]))/length(genes[[i]])
    }
  }
  cl = hclust(dist(resmat), method = "ward.D2")
  cl = cutree(cl, ncl)
  res_df = data.frame("Description" = df$Description, "qvalue" = df$qvalue, cl, 
                      "geneID" = unlist(df$geneID), stringsAsFactors = F)
  res_df = res_df[order(res_df$qvalue, decreasing = F),]
  topterms = unlist(tapply(res_df$Description, res_df$cl, function(x) x[1:nt]))
  res_df = res_df[res_df$Description %in% topterms,]
  
  return(res_df)
}

getTopTermsPaired = function(godf, genescol = "genes_all", ncl = 5, nt = 2){
  if(nrow(godf)<ncl) return(godf)
  genes = sapply(godf[,genescol], function(x) strsplit(x, "/"))
  resmat = matrix(0, length(genes), length(genes))
  for(i in 1:length(genes)){
    for(j in 1:length(genes)){
      resmat[i,j] = length(intersect(genes[[i]], genes[[j]]))/length(genes[[i]])
    }
  }
  cl = hclust(dist(resmat), method = "ward.D2")
  cl = cutree(cl, ncl)
  res_df = data.frame("Description" = godf$Description, "qval_embolized" = godf$qval_embolized, 
                      "qval_regenerating" = godf$qval_regenerating, 
                      "genes_embolized" = unlist(godf$genes_embolized), 
                      "genes_regenerating" = unlist(godf$genes_regenerating),
                      cl, stringsAsFactors = F)
  res_df$qval_max = apply(res_df[,c("qval_embolized","qval_regenerating")], 1, function(x) min(x))
  res_df = res_df[order(res_df$qval_max, decreasing = F),]
  topterms = unlist(tapply(res_df$Description, res_df$cl, function(x) x[1:nt]))
  res_df = res_df[res_df$Description %in% topterms,]
  
  return(res_df)
}

rpfilter = function(x){
  return(!grepl("^RPL", x$geneID) & !grepl("^RPS", x$geneID) & 
           !grepl("/RPL", x$geneID, fixed = T) & !grepl("/RPS", x$geneID, fixed = T))
}

Figure 1

Main Figure

UMAP with all cell types

colsmajor = c("Cholangiocytes" = "bisque2", "Endothelial" = "aquamarine4",  
              "Hepatocytes" = "tomato3", "Immune" = "skyblue", 
              "Mesenchymal" = "sandybrown", "Doublets" = "grey90")
colcond = c("healthy" = "orange", "regenerating" = "salmon", 
            "embolised" = "darkred", "embolized" = "darkred")
coldon = c("sc_H1" = "plum4", "sc_H2" = "salmon4", "sc_H3" = "lightsalmon3")
coldonall = c("sc_H1" = "plum4", "sc_H2" = "salmon4", "sc_H3" = "lightsalmon3", 
              "sc_E1" = "darkorange3", "sc_R1" = "darkorange3", "sc_E2" = "goldenrod3", 
              "sc_R2" = "goldenrod3",  "sc_E3" = "mediumorchid", "sc_R3" = "mediumorchid", 
              "sc_E4" = "plum2", "sc_R4" = "plum2", "sc_E5" = "palevioletred3", 
              "sc_R5" = "palevioletred3", "sc_E6" = "peachpuff2", "sc_R6" = "peachpuff2",
              "sc_E1/sc_R1" = "darkorange3", "sc_E2/sc_R2" = "goldenrod3", 
              "sc_E3/sc_R3" = "mediumorchid", "sc_E4/sc_R4" = "plum2", 
              "sc_E5/sc_R5" = "palevioletred3", "sc_E6/sc_R6" = "peachpuff2")

colsallct = c("Cholangiocytes" = "bisque2", "Hepatocytes" = "tomato3", "Stellate cells" = "sandybrown", 
              "Doublets" = "grey90", "LSEC pericentral" = "aquamarine4", "LSEC periportal" = "aquamarine2", 
              "Endothelial cells (non-LSEC)" = "forestgreen", "Plasmablasts" = "darkorchid4",
              "Kupffer cells" = "skyblue1", "ab-T cells" = "lightskyblue3", "gd-T cells" = "cornflowerblue",
              "B cells" = "darkorchid1", "cDCs" = "deepskyblue1", "pDCs" = "royalblue3", 
              "Macrophages" = "steelblue3", "Dividing cells" = "grey35", "LSEC" = "aquamarine3")
colsub = c("rosybrown4","thistle4","thistle3","aquamarine4","aquamarine3",
           "forestgreen","tomato4","tomato3","darksalmon","skyblue",
           "cadetblue","sandybrown","palegoldenrod")
colsmidct = c("Cholangiocytes" = "bisque2", "Hepatocytes" = "tomato3", 
              "Mesenchymal" = "sandybrown", "pDCs" = "royalblue3", 
              "T cells" = "lightskyblue3", "LSEC" = "aquamarine3",
              "other ECs" = "forestgreen", "B cells" = "darkorchid1",
              "Kupffer cells" = "skyblue1", "other Mono-Mac" = "steelblue3", 
              "ILC" = "slateblue1")

plot_df = data.frame(hcells_css@reductions$umap_css@cell.embeddings)
plot_df$names_major = as.character(hcells_css@meta.data$names_major)
plot_df$names_major[plot_df$names_major=="T cells" | 
                      plot_df$names_major=="B cells" | 
                      plot_df$names_major=="Macrophages" |
                      plot_df$names_major=="DCs"] = "Immune"
plot_df$names_major[plot_df$names_major=="Stellate cells"] = "Mesenchymal"
plot_df$names_major[plot_df$names_major=="Endothelial cells"] = "Endothelial"
plt = ggplot(plot_df, aes(x = UMACSS_1, y = UMACSS_2, colour = names_major))+
  geom_point(size = pointsize)+
  guides(colour = guide_legend(override.aes = list(size = 3), title = "Cell Type"))+
  scale_colour_manual(values = colsmajor)+
  theme_classic()+
  th_gen+
  theme(axis.line = element_blank(),
        axis.ticks = element_blank(),
        axis.title = element_blank(),
        axis.text = element_blank(),
        legend.title = element_text(hjust = 0),
        aspect.ratio = 1)

pdf("figure_panels/fig1/umap_fresh_major.pdf", 
    useDingbats = F, height = 4, width = 5)
print(plt)
dev.off()
png("figure_panels/fig1/umap_fresh_major.png", 
    height=8, width=8, unit="cm", res=600, antialias = "subpixel")
print(plt)
dev.off()
pdf("figure_panels/fig1/umap_fresh_major_noLeg.pdf", 
    useDingbats = F, height = 4, width = 5)
print(plt+theme(legend.position = "none"))
dev.off()
png("figure_panels/fig1/umap_fresh_major_noLeg.png", 
    height=8, width=8, unit="cm", res=600, antialias = "subpixel")
print(plt+theme(legend.position = "none"))

UMAP with healthy donors

plot_df = data.frame(hcells_css@reductions$umap_css@cell.embeddings)
plot_df$Donor = factor(hcells_css@meta.data$Donor)
plot_df$Donor = plyr::revalue(plot_df$Donor, c("HD1" = "sc_H1", "HD2" = "sc_H2", "HD3" = "sc_H3"))

plt = ggplot(plot_df[sample(1:nrow(plot_df), nrow(plot_df), replace = F),], 
             aes(x = UMACSS_1, y = UMACSS_2, colour = Donor))+
  geom_point(size = pointsize)+
  guides(colour = guide_legend(override.aes = list(size = 3), title = "Donors"))+
  scale_colour_manual(values = coldon)+
  theme_classic()+
  th_gen+
  theme(axis.line = element_blank(),
        axis.ticks = element_blank(),
        axis.title = element_blank(),
        axis.text = element_blank(),
        legend.title = element_text(hjust = 0),
        aspect.ratio = 1)

pdf("figure_panels/fig1/umap_fresh_donors.pdf", 
    useDingbats = F, height = 4, width = 5)
print(plt)
dev.off()
png("figure_panels/fig1/umap_fresh_donors.png", 
    height=8, width=8, unit="cm", res=600, antialias = "subpixel")
print(plt)
dev.off()
pdf("figure_panels/fig1/umap_fresh_donors_noLeg.pdf", 
    useDingbats = F, height = 4, width = 5)
print(plt+theme(legend.position = "none"))
dev.off()
png("figure_panels/fig1/umap_fresh_donors_noLeg.png", 
    height=8, width=8, unit="cm", res=600, antialias = "subpixel")
print(plt+theme(legend.position = "none"))
dev.off()

Violins for cell type markers

markers = c("ASGR1", "APOB", "APOC3",
            "KRT7","CLDN4","EPCAM",
            "CLEC4G","PECAM1","CD36",
            "DCN","COLEC11","ACTA2",
            "PTPRC","HLA-DQA1","LYZ")

exp_mk = data.frame(Matrix::t(hcells_css@assays$SCT@data[markers,]))
plot_df = cbind(exp_mk, 
                data.frame(names_major = as.character(hcells_css@meta.data[,c("names_major")]), 
                           stringsAsFactors = F))
plot_df$names_major[plot_df$names_major=="T cells" | 
                      plot_df$names_major=="B cells" | 
                      plot_df$names_major=="Macrophages" | 
                      plot_df$names_major=="DCs"] = "Immune"
plot_df$names_major[plot_df$names_major=="Endothelial cells"] = "Endothelial"
plot_df$names_major[plot_df$names_major=="Stellate cells"] = "Mesenchymal"
plot_df = plot_df[plot_df$names_major!="Doublets",]
plot_df$names_major = factor(plot_df$names_major, 
                             levels = c("Hepatocytes", "Cholangiocytes", "Endothelial", "Mesenchymal",
                                        "Immune"))
plot_df = reshape2::melt(plot_df)
plot_df$variable = factor(gsub(".", "-", plot_df$variable, fixed = T), levels = markers)

vio_mk = ggplot(plot_df, aes(x = names_major, y = value, fill = names_major))+
  facet_grid(variable~names_major, scales = "free")+
  geom_violin(scale = "width", size = 0.3)+
  scale_y_continuous(breaks = seq(0,10,2), labels = seq(0,10,2), name = "log(exp+1)")+
  scale_fill_manual(values = colsmajor)+
  theme(strip.background.y = element_blank(),
        strip.background.x = element_rect(fill = "transparent", 
                                          colour = "black", size = 0.8),
        strip.text.y = element_text(angle = 0, size = 6.5, colour = "black", face = "bold"),
        strip.text.x = element_text(face = "bold", size = 6.5, colour = "black"),
        legend.position = "none",
        axis.text.x = element_blank(),
        axis.title.x = element_blank(), 
        axis.line.x.bottom = element_blank())

pdf("figure_panels/fig1/violin_markers_major.pdf", 
    useDingbats = F, height = 5, width = 5.5)
print(vio_mk)
dev.off()
png("figure_panels/fig1/violin_markers_major.png", height = 425, width = 450, antialias = "subpixel")
print(vio_mk)
dev.off()

Heatmap for major cell types

# removed HHIP
nmg=c("PDGFRA","CALD1","COL6A1","PDGFRB",
      "CSF1R","CD163","MARCO","CD69","IL7R","PCK1","CYP2A7","CYP3A4","CRP",
      "ROBO4","EGFL7","CLEC4M","FCN2","KRT7","CFTR","ONECUT1")
plot_df = data.frame(row.names = rownames(hcells_css@meta.data),
                     names_major = as.character(hcells_css@meta.data[,c("names_major")]),
                     donor = as.character(hcells_css@meta.data[,c("Donor")]),
                     stringsAsFactors = F) 
plot_df$names_major[plot_df$names_major=="T cells" | 
                      plot_df$names_major=="B cells" | 
                      plot_df$names_major=="Macrophages" | 
                      plot_df$names_major=="DCs"] = "Immune"
plot_df$names_major[plot_df$names_major=="Endothelial cells"] = "Endothelial"
plot_df$names_major[plot_df$names_major=="Stellate cells"] = "Mesenchymal"
plot_df = plot_df[plot_df$names_major!="Doublets",]
plot_df = plot_df[order(plot_df$names_major),]
m1=GetAssayData(hcells_css, slot="data")[nmg,rownames(plot_df)]
m1 = t(apply(m1, 1, scale, scale = F))
colnames(m1)=rep("",ncol(m1))
coul <- colorRampPalette(brewer.pal(9, "Greys"))(100)[-c(1:5)]
m1[m1>1] = 1
#m1[m1<=(-2)] = -2
heatmap(m1,Rowv = NA,Colv = NA, col=coul)

Violins for all cell types

exp_mk = data.frame(Matrix::t(hcells_css@assays$SCT@data[markers,]))
plot_df = cbind(exp_mk, 
                hcells_css@meta.data[,c("names_major","names_clusters")])
plot_df$names_major = as.character(plot_df$names_major)
plot_df$names_major[plot_df$names_major=="T cells" | 
                      plot_df$names_major=="B cells" | 
                      plot_df$names_major=="Macrophages" | 
                      plot_df$names_major=="DCs"] = "Immune"
plot_df$names_major[plot_df$names_major=="Endothelial cells"] = "Endothelial"
plot_df$names_major[plot_df$names_major=="Stellate cells"] = "Mesenchymal"
plot_df = plot_df[plot_df$names_major!="Doublets",]
plot_df$names_major = factor(plot_df$names_major, 
                             levels = c("Hepatocytes", "Cholangiocytes", "Endothelial", "Mesenchymal",
                                        "Immune"))
plot_df = reshape2::melt(plot_df)
plot_df$variable = factor(gsub(".", "-", plot_df$variable, fixed = T), levels = markers)

ggplot(plot_df, aes(x = names_clusters, y = value, fill = names_clusters))+
  facet_grid(variable~names_major, scales = "free")+
  geom_violin(scale = "width")+
  scale_y_continuous(breaks = seq(0,10,2), labels = seq(0,10,2), name = "log(exp+1)")+
  #scale_fill_manual(values = colsmajor)+
  theme(strip.background.y = element_blank(),
        strip.background.x = element_rect(fill = "transparent", 
                                          colour = "black", size = 0.8),
        strip.text.y = element_text(angle = 0, size = 6.5, colour = "black", face = "bold"),
        strip.text.x = element_text(face = "bold", size = 6.5, colour = "black"),
        legend.position = "none",
        axis.text.x = element_text(angle = -35, hjust = 0, vjust = 0.1),
        axis.line.x.bottom = element_blank())

Cell type proportions

plot_df = data.frame("names_major" = as.character(hcells_css@meta.data$names_major),
                     "Donor" = as.character(hcells_css@meta.data$Donor),
                     stringsAsFactors = F)
plot_df$names_major[plot_df$names_major=="T cells" | 
                      plot_df$names_major=="B cells" | 
                      plot_df$names_major=="Macrophages" | 
                      plot_df$names_major=="DCs"] = "Immune"
plot_df$names_major[plot_df$names_major=="Endothelial cells"] = "Endothelial"
plot_df$names_major[plot_df$names_major=="Stellate cells"] = "Mesenchymal"
plot_df = plot_df[plot_df$names_major!="Doublets",]
cnts_ct = table(plot_df$names_major, plot_df$Donor)
plot_df = reshape2::melt(apply(cnts_ct, 1, function(x) x/sum(x)))
#plot_df$Var2 = factor(plot_df$Var2, 
#                      levels = c("Hepatocytes", "Cholangiocytes", "Endothelial","Mesenchymal", "Immune"))
plot_df$Var1 = plyr::revalue(plot_df$Var1, c("HD1" = "sc_H1", "HD2" = "sc_H2", "HD3" = "sc_H3"))
plot_df$Var1 = factor(plot_df$Var1, levels = rev(levels(plot_df$Var1)))

plt = ggplot(plot_df, aes(x = Var2, y = value*100, fill = Var1))+
  geom_bar(stat = "identity")+
  scale_y_continuous(expand = c(0,0))+
  scale_fill_manual(values = coldon)+
  labs(y = "Cell type proportion [%]", x = NULL)+
  theme_classic()+
  th_gen+
  theme(axis.line.x = element_blank(),
        axis.ticks.y = element_line(),
        axis.ticks.x = element_blank())

pdf("figure_panels/fig1/proportions_fresh_major_donor.pdf", 
    useDingbats = F, height = 4, width = 5)
print(plt)
dev.off()
png("figure_panels/fig1/proportions_fresh_major_donor.png", 
    height = 325, width = 400, antialias = "subpixel")
print(plt)
dev.off()
pdf("figure_panels/fig1/proportions_fresh_major_donor_noLeg.pdf", 
    useDingbats = F, height = 4, width = 5)
print(plt+theme(legend.position = "none"))
dev.off()
png("figure_panels/fig1/proportions_fresh_major_donor_noLeg.png", 
    height = 325, width = 400, antialias = "subpixel")
print(plt+theme(legend.position = "none"))
dev.off()

Import markers

markers = c("ASGR1", "APOB", "APOC3", "CYP2E1", "HAMP", "ORM1", "SAA1", "NNMT", "FABP1", "MT1G", "ORM2", "TTR", "HP", "APOC1", "APOA2","FGB","CYP3A4","CPS1","ARG1","SAA2",# Hep
            "KRT7","CLDN4","EPCAM", "TACSTD2", "CD24", "KRT19", "ANXA4", "CXCL6", "FXYD2", "SOX4",  "CRYAB","DEFB1","SLC12A2","MMP7","TNFRSF12A","CXCL1","BICC1","S100A14","DCDC2","PLPP2",# Cho
            "CLEC4G","PECAM1","CD36", "FCN3", "DNASE1L3", "CLEC1B", "CRHBP", "AKAP12", "IFI27", "GNG11", "IL33","FLT1","PRSS23","ENG","RAMP3","F8","VWF","CLDN5","CCL14","LYVE1",# End
            "DCN","COLEC11","ACTA2", "CCL2", "TAGLN", "IGFBP3", "BGN", "LUM", "COL3A1", "MYL9", "TPM2","AEBP1","GGT5","ASPN","COL14A1","PTGDS","COL6A1","CYR61","COLEC10","CXCL12",# Mes
            "PTPRC","HLA-DQA1","LYZ", "LILRB2", "MARCO", "C1QB", "FCGR3A",  "NKG7", "MS4A1", "MZB1", "CCL5","KLRD1","AREG","MS4A7","AXL","CD69","GPR183","TLR2","CD44","IL7R") # Imm
mk_list = list("Hepatocytes" = markers[1:10], "Cholangiocytes" = markers[21:30], 
               "Endothelial" = markers[41:50], "Mesenchymal" = markers[61:70], "Immune" = markers[81:90])

cell_type_mk = readRDS(file = "results/integr_allcells/cell_type_mk_major.RDS")
fresh_de = list()
for(n in unique(cell_type_mk$major_all$cluster)){
  fresh_de[[n]] = cell_type_mk$major_all[cell_type_mk$major_all$cluster==n,]
}

load("data/processed/received_Aga/DE_tables_celltype_frozenSCT.rdata")
frozen_de = list("Cholangiocytes" = cell_type_mk_cho, "Hepatocytes" = cell_type_mk_hep, 
                 "Endothelial" = cell_type_mk_ec, "Immune" = cell_type_mk_imm, 
                 "Mesenchymal" = cell_type_mk_mes)
for(n in names(frozen_de)){
  frozen_de[[n]]$cluster = n
  frozen_de[[n]]$gene = rownames(frozen_de[[n]])
}

both_list = list()
top_mk = list()
for(n in names(fresh_de)){
  both_DE = merge(fresh_de[[n]], frozen_de[[n]], by = "gene", all = T)[,c(1,3,6,9,12)]
  colnames(both_DE) = c("gene", "FC_fresh", "pval_fresh", "FC_frozen", "pval_frozen")
  both_DE$FC_fresh[is.na(both_DE$FC_fresh)] = 0
  both_DE$pval_fresh[is.na(both_DE$pval_fresh)] = 1
  both_DE$FC_frozen[is.na(both_DE$FC_frozen)] = 0
  both_DE$pval_frozen[is.na(both_DE$pval_frozen)] = 1
  both_DE$celltype = n
  both_DE$s = both_DE$FC_fresh+both_DE$FC_frozen
  both_DE = both_DE[order(both_DE$s, decreasing = T),]
  cond_fc = both_DE$FC_fresh>=0.1 & both_DE$FC_frozen>=0.1
  cond_pv = both_DE$pval_fresh<=0.05 | both_DE$pval_frozen<=0.05
  both_DE$iscol = ifelse(cond_fc & cond_pv, 
                         ifelse(both_DE$gene %in% mk_list[[n]], "istop", "isDE"), "notDE")
  both_DE$iscol = factor(both_DE$iscol, levels = rev(c("istop", "isDE", "notDE")))
  both_DE$istop = both_DE$gene %in% mk_list[[n]]#[1:3]
  both_list[[n]] = both_DE
  top_mk[[n]] = both_DE$gene[cond_fc & cond_pv][1:20]
}
top_mk = top_mk[c(4,1,3,5,2)]

Heatmap per cluster

cell_type_mk = readRDS(file = "results/cond_effect/cell_type_mk.RDS")
# cut markers here depending on the size you want for heatmap, then adjust figure
markers = c("ASGR1", "APOB", "APOC3", "CYP2E1", "HAMP", "ORM1", "SAA1", "NNMT", "FABP1", "MT1G", "ORM2", "TTR", "HP", "APOC1", "APOA2","FGB","CYP3A4","CPS1","ARG1","SAA2",# Hep
            "KRT7","CLDN4","EPCAM", "TACSTD2", "CD24", "KRT19", "ANXA4", "CXCL6", "FXYD2", "SOX4",  "CRYAB","DEFB1","SLC12A2","MMP7","TNFRSF12A","CXCL1","BICC1","S100A14","DCDC2","PLPP2",# Cho
            "CLEC4G","PECAM1","CD36", "FCN3", "DNASE1L3", "CLEC1B", "CRHBP", "AKAP12", "IFI27", "GNG11", "IL33","FLT1","PRSS23","ENG","RAMP3","F8","VWF","CLDN5","CCL14","LYVE1",# End
            "DCN","COLEC11","ACTA2", "CCL2", "TAGLN", "IGFBP3", "BGN", "LUM", "COL3A1", "MYL9", "TPM2","AEBP1","GGT5","ASPN","COL14A1","PTGDS","COL6A1","CYR61","COLEC10","CXCL12",# Mes
            "PTPRC","HLA-DQA1","LYZ", "LILRB2", "MARCO", "C1QB", "FCGR3A",  "NKG7", "MS4A1", "MZB1", "CCL5","KLRD1","AREG","MS4A7","AXL","CD69","GPR183","TLR2","CD44","IL7R") # Imm
markers_topboth = unlist(top_mk)
markers = markers_topboth

load("data/processed/received_Aga/mean_exp_celltype_frozen.rdata")
hcells_css@meta.data$major_ct = as.character(hcells_css@meta.data$names_major)
hcells_css@meta.data$major_ct[hcells_css@meta.data$major_ct=="T cells" | 
                      hcells_css@meta.data$major_ct=="B cells" | 
                      hcells_css@meta.data$major_ct=="Macrophages" |
                      hcells_css@meta.data$major_ct=="DCs"] = "Immune"
hcells_css@meta.data$major_ct[hcells_css@meta.data$major_ct=="Stellate cells"] = "Mesenchymal"
hcells_css@meta.data$major_ct[hcells_css@meta.data$major_ct=="Endothelial cells"] = "Endothelial"
m_ct_fresh = apply(as.matrix(hcells_css@assays$SCT@data), 1, 
                   function(x) tapply(x, hcells_css@meta.data$major_ct, mean))

m_ct_frozen = t(m_ct_frozen)
m_ct_frozen = t(apply(m_ct_frozen, 1, scale))
colnames(m_ct_frozen) = rownames(m_ct_fresh)[-2]
m_ct_fresh = t(m_ct_fresh)[,-2]
m_ct_fresh = t(apply(m_ct_fresh, 1, scale))
colnames(m_ct_fresh) = colnames(m_ct_frozen)

mean_ct = merge(m_ct_fresh, m_ct_frozen, by = 0)
rownames(mean_ct) = mean_ct[,1]
mean_ct = mean_ct[,-1]
colnames(mean_ct) = paste0(rep(colnames(m_ct_fresh), 2), rep(c("_fresh", "_frozen"), each = 5))

hcl = hclust(dist(t(mean_ct)), method = "ward.D2")
plot(hcl)

annot_df = data.frame(row.names = colnames(mean_ct),
                      "cell type" = unlist(lapply(strsplit(colnames(mean_ct), "_"), function(x) x[1])),
                      "processing" = unlist(lapply(strsplit(colnames(mean_ct), "_"), function(x) x[2])))

callback = function(hc, mat){
  sv = svd(t(mat))$v[,2]
    dend = reorder(as.dendrogram(hc), c(1:4,))
    as.hclust(dend)
}

pdf("figure_panels/fig1/major_celltypes_heatmap_freshfrozen_10.pdf", 
    useDingbats = F, height = 7.5, width = 6)
heat = pheatmap::pheatmap(mean_ct[markers[c(41:50, 61:70, 81:90,  1:10, 21:30)],], 
                          clustering_method = "ward.D2", treeheight_row = 0, annotation_col = annot_df,
                          cluster_rows = F, treeheight_col = 20, fontsize_row = 5.5,
                          annotation_colors = list("processing" = c("fresh" = "#d4d4d4", 
                                                                    "frozen" = "#5d5d5d"),
                                                   "cell.type" = colsmajor[-6]),
                          show_colnames = F, 
                          color = colorRampPalette(rev(brewer.pal(n = 9, name = "RdBu")))(100))
dev.off()

FC fresh vs frozen

mk_list = list("Hepatocytes" = markers[1:10], "Cholangiocytes" = markers[21:30], 
               "Endothelial" = markers[41:50], "Mesenchymal" = markers[61:70], "Immune" = markers[81:90])

plt_fc_list = list()
pdf("figure_panels/fig1/major_celltypes_FCscatter_10.pdf", 
    useDingbats = F, height = 5, width = 4.5)
for(n in names(both_list)){
  both_DE = both_list[[n]]
  cc = cor(both_DE$FC_fresh, both_DE$FC_frozen)
  both_DE = both_DE[order(both_DE$iscol),]
  
  cols_use = c("istop" = unname(colsmajor[n]), "isDE" = "grey35", "notDE" = "grey85")
  plt = ggplot(both_DE, aes(x = FC_fresh, y = FC_frozen, colour = iscol))+
    geom_vline(xintercept = 0, colour = "grey40")+
    geom_hline(yintercept = 0, colour = "grey40")+
    geom_point()+
    geom_text_repel(data = both_DE[both_DE$istop & both_DE$iscol=="istop",], mapping = aes(label = gene),
                    fontface = "bold")+
    theme_bw()+
    scale_colour_manual(values = cols_use)+
    labs(title = n, subtitle = paste0("PCC = ", round(cc, 2)))+
    theme(aspect.ratio = 1,
          axis.text = element_text(colour = "black"),
          legend.position = "none", 
          panel.grid = element_blank())
  
  plt_fc_list[[n]] = plt
  print(plt)
}
dev.off()

Expression comparison fresh vs frozen

load("data/processed/received_Aga/mean_exp_celltype_frozenCOUNTS.rdata")
hcells_css@meta.data$major_ct = as.character(hcells_css@meta.data$names_major)
hcells_css@meta.data$major_ct[hcells_css@meta.data$major_ct=="T cells" | 
                      hcells_css@meta.data$major_ct=="B cells" | 
                      hcells_css@meta.data$major_ct=="Macrophages" |
                      hcells_css@meta.data$major_ct=="DCs"] = "Immune"
hcells_css@meta.data$major_ct[hcells_css@meta.data$major_ct=="Stellate cells"] = "Mesenchymal"
hcells_css@meta.data$major_ct[hcells_css@meta.data$major_ct=="Endothelial cells"] = "Endothelial"
m_ct_fresh = apply(as.matrix(hcells_css@assays$SCT@counts), 1, 
                   function(x) tapply(x, hcells_css@meta.data$major_ct, mean))

m_ct_frozen = log(t(m_ct_frozen.counts))
#m_ct_frozen = t(apply(m_ct_frozen, 1, scale))
colnames(m_ct_frozen) = rownames(m_ct_fresh)[-2]
m_ct_fresh = log(t(m_ct_fresh)[,-2])
#m_ct_fresh = t(apply(m_ct_fresh, 1, scale))
colnames(m_ct_fresh) = colnames(m_ct_frozen)

mean_ct = merge(m_ct_fresh, m_ct_frozen, by = 0)
rownames(mean_ct) = mean_ct[,1]
mean_ct = mean_ct[,-1]
colnames(mean_ct) = paste0(rep(colnames(m_ct_fresh), 2), rep(c("_fresh", "_frozen"), each = 5))

par(mfrow = c(2,3))
plot(mean_ct$Cholangiocytes_fresh, mean_ct$Cholangiocytes_frozen, pch = 20, cex = 0.5, xlim = c(-8,5), ylim = c(-8,5))
title("Cholangiocytes")
abline(0,1)
plot(mean_ct$Endothelial_fresh, mean_ct$Endothelial_frozen, pch = 20, cex = 0.5, xlim = c(-8,5), ylim = c(-8,5))
title("Endothelial")
abline(0,1)
plot(mean_ct$Hepatocytes_fresh, mean_ct$Hepatocytes_frozen, pch = 20, cex = 0.5, xlim = c(-10,6), ylim = c(-10,6))
title("Hepatocytes")
abline(0,1)
plot(mean_ct$Immune_fresh, mean_ct$Immune_frozen, pch = 20, cex = 0.5, xlim = c(-8,5), ylim = c(-8,5))
title("Immune")
abline(0,1)
plot(mean_ct$Mesenchymal_fresh, mean_ct$Mesenchymal_frozen, pch = 20, cex = 0.5, xlim = c(-8,5), ylim = c(-8,5))
title("Mesenchymal")
abline(0,1)
plot(rowMeans(mean_ct[,1:5]), rowMeans(mean_ct[,6:10]), pch = 20, cex = 0.5, xlim = c(-8,5), ylim = c(-8,5))
title("All")
abline(0,1)

Supplementary

UMAP with healthy donors split

plot_df = data.frame(hcells_css@reductions$umap_css@cell.embeddings)
plot_df$Donor = factor(hcells_css@meta.data$Donor)
plot_df$Donor = plyr::revalue(plot_df$Donor, c("HD1" = "sc_H1", "HD2" = "sc_H2", "HD3" = "sc_H3"))

plt = ggplot(plot_df[sample(1:nrow(plot_df), nrow(plot_df), replace = F),], 
             aes(x = UMACSS_1, y = UMACSS_2, colour = Donor))+
  facet_wrap(~Donor)+
  geom_point(size = 0.15)+
  guides(colour = guide_legend(override.aes = list(size = 3), title = "Donors"))+
  scale_colour_manual(values = coldon)+
  theme_classic()+
  th_gen+
  theme(axis.line = element_blank(),
        axis.ticks = element_blank(),
        axis.title = element_blank(),
        axis.text = element_blank(),
        legend.title = element_text(hjust = 0),
        aspect.ratio = 1)

pdf("figure_panels/fig1/umap_fresh_donors_split.pdf", 
    useDingbats = F, height = 4, width = 7)
print(plt)
dev.off()
png("figure_panels/fig1/umap_fresh_donors_split.png", 
    height=8, width=10, unit="cm", res=600, antialias = "subpixel")
print(plt)
dev.off()
pdf("figure_panels/fig1/umap_fresh_donors_noLeg_split.pdf", 
    useDingbats = F, height = 4, width = 7)
print(plt+theme(legend.position = "none"))
dev.off()
png("figure_panels/fig1/umap_fresh_donors_noLeg_split.png", 
    height=8, width=10, unit="cm", res=600, antialias = "subpixel")
print(plt+theme(legend.position = "none"))
dev.off()

Boxplot with nUMI

plot_df = data.frame(hcells_css@reductions$umap_css@cell.embeddings)
plot_df$names_major = as.character(hcells_css@meta.data$names_major)
plot_df$names_major[plot_df$names_major=="T cells" | 
                      plot_df$names_major=="B cells" | 
                      plot_df$names_major=="Macrophages" |
                      plot_df$names_major=="DCs"] = "Immune"
plot_df$names_major[plot_df$names_major=="Stellate cells"] = "Mesenchymal"
plot_df$names_major[plot_df$names_major=="Endothelial cells"] = "Endothelial"
plot_df$nGene = hcells_css@meta.data$nFeature_SCT
plot_df = plot_df[,3:5]
saveRDS(plot_df, "./to_send/df_celltypes_umi_gene.RDS") # plotted by Aga

UMAP with all cell types

plot_df = data.frame(hcells_css@reductions$umap_css@cell.embeddings)
colnames(plot_df) = c("UMAPCSS_1", "UMAPCSS_2")
plot_df$names_major = as.character(hcells_css@meta.data$names_major)
plot_df$names_major[plot_df$names_major=="T cells" | 
                      plot_df$names_major=="B cells" | 
                      plot_df$names_major=="Macrophages" |
                      plot_df$names_major=="DCs"] = "Immune"
plot_df$names_major[plot_df$names_major=="Stellate cells"] = "Mesenchymal"
plot_df$names_major[plot_df$names_major=="Endothelial cells"] = "Endothelial"
plot_df$names_clusters = as.character(hcells_css@meta.data$names_clusters)
plot_df$names_clusters[grepl("Hepatocytes ", plot_df$names_clusters)] = "Hepatocytes"
plot_df$names_clusters[grepl("central vein", plot_df$names_clusters)] = "LSEC pericentral"
plot_df$names_clusters[grepl("portal vein", plot_df$names_clusters)] = "LSEC periportal"
plot_df$names_clusters[grepl("interaction", plot_df$names_clusters) |
                         grepl("mix", plot_df$names_clusters) ] = "Doublets"
plot_df$names_clusters = factor(plot_df$names_clusters, 
                                levels = c("Cholangiocytes", "LSEC periportal", "LSEC pericentral", 
                                           "Endothelial cells (non-LSEC)", "Hepatocytes", "Kupffer cells",
                                           "Macrophages", "cDCs", "pDCs", "ab-T cells", "gd-T cells",
                                           "B cells", "Plasmablasts",
                                           "Stellate cells", "Dividing cells", "Doublets"))

pltumap = ggplot(plot_df, aes(y = UMAPCSS_2, x = UMAPCSS_1, colour = names_clusters))+
  geom_point(size = pointsize)+
  labs(colour = "Cell Type")+
  scale_colour_manual(values = c(colsallct, "Other" = "grey88"))+
  guides(colour = guide_legend(override.aes = list(size = 2)))+
  theme(aspect.ratio = 1,
        legend.title = element_text(hjust = 0),
        axis.line = element_blank(),
        axis.text = element_blank())

pdf("figure_panels/fig1/umap_h_all_celltypes.pdf", height = 5, width = 6, useDingbats = F)
print(pltumap)
dev.off()
png("figure_panels/fig1/umap_h_all_celltypes.png", height = 325, width = 400, antialias = "subpixel")
print(pltumap)
dev.off()

plot_df$names_clusters[plot_df$names_major!="Immune"] = "Doublets"
levels(plot_df$names_clusters)[levels(plot_df$names_clusters)=="Doublets"] = "Other"
pltumap = ggplot(plot_df, aes(y = UMAPCSS_2, x = UMAPCSS_1, colour = names_clusters))+
  geom_point(size = pointsize)+
  labs(colour = "Cell Type")+
  scale_colour_manual(values = c(colsallct, "Other" = "grey88"))+
  guides(colour = guide_legend(override.aes = list(size = 2)))+
  theme(aspect.ratio = 1,
        legend.title = element_text(hjust = 0),
        axis.line = element_blank(),
        axis.text = element_blank())

Heatmap with markers for all cell types

plot_df = data.frame(hcells_css@reductions$umap_css@cell.embeddings)
colnames(plot_df) = c("UMAPCSS_1", "UMAPCSS_2")
plot_df$names_major = as.character(hcells_css@meta.data$names_major)
plot_df$names_major[plot_df$names_major=="T cells" | 
                      plot_df$names_major=="B cells" | 
                      plot_df$names_major=="Macrophages" |
                      plot_df$names_major=="DCs"] = "Immune"
plot_df$names_major[plot_df$names_major=="Stellate cells"] = "Mesenchymal"
plot_df$names_major[plot_df$names_major=="Endothelial cells"] = "Endothelial"
plot_df$names_clusters = as.character(hcells_css@meta.data$names_clusters)
plot_df$names_clusters[grepl("Hepatocytes ", plot_df$names_clusters)] = "Hepatocytes"
plot_df$names_clusters[grepl("central vein", plot_df$names_clusters)] = "LSEC pericentral"
plot_df$names_clusters[grepl("portal vein", plot_df$names_clusters)] = "LSEC periportal"
plot_df$names_clusters[grepl("interaction", plot_df$names_clusters) |
                         grepl("mix", plot_df$names_clusters) ] = "Doublets"
plot_df$names_clusters = factor(plot_df$names_clusters, 
                                levels = c("Cholangiocytes", "LSEC periportal", "LSEC pericentral", 
                                           "Endothelial cells (non-LSEC)", "Hepatocytes", "Kupffer cells",
                                           "Macrophages", "cDCs", "pDCs", "ab-T cells", "gd-T cells",
                                           "B cells", "Plasmablasts",
                                           "Stellate cells", "Dividing cells", "Doublets"))
plot_df = plot_df[plot_df$names_clusters!="Doublets",]

nmg = c("KRT7", "CD24",  "LILRA4", "IRF7",  "FCER1A", "CLEC10A",  "CCL23","CLEC1B",  "INMT","SELP",  
        "MARCO","CD5L",  "CYP2E1","ORM2",  "MS4A1","CD79A",  "CLEC14A","EDNRB",  "IGLC2","CD27",  
        "TRDC","KLRD1",  "TRAC","CD3E",  "DCN","COLEC11",  "S100A12","CD300E")
nmg = nmg[c(1,2,13,14,17,18,7,8,9,10,25,26,27,28,11,12,5,6,3,4,23,24,21,22,15,16,19,20)]

m1 = GetAssayData(hcells_css, slot="data")[nmg,rownames(plot_df)]
m_ct_fresh = apply(m1, 1, function(x) tapply(x, as.character(plot_df$names_clusters), mean))
m_ct_fresh[m_ct_fresh>4] = 4.01
m_ct_fresh[m_ct_fresh<(-4)] = -4.01

annot_df = unique(plot_df[,c(3,4)])
annot_df = data.frame(row.names = annot_df$names_clusters, "ct" = annot_df$names_major)
heat = pheatmap::pheatmap(scale(m_ct_fresh[c(4,7,10,9,5,14,11,8,3,12,1,6,2,13),]), 
                          clustering_method = "ward.D2", treeheight_row = 0,
                          treeheight_col = 20, fontsize_row = 8, fontsize_col = 8,
                          cluster_cols = F, cluster_rows = F,
                          annotation_row = annot_df, annotation_colors = list(ct = colsmajor[-6]),
                          color = colorRampPalette(rev(brewer.pal(n = 9, name = "RdBu")))(100))

pdf("figure_panels/fig1/heatmap_all_celltypes.pdf", height = 3.5, width = 7, useDingbats = F)
print(heat)
dev.off()

Figure 2

Main Figure

UMAP with cell types

allcells_css = readRDS(file = "data/processed/allcells_css_reannot.RDS")

load(file = "data/processed/healthy_srat_beforeFiltering.RData")
load(file = "data/processed/cond_srat_beforeFiltering.RData")

allcells_css = readRDS(file = "data/processed/allcells_css_reannot.RDS")

load(file = "data/processed/healthy_srat_beforeFiltering.RData")
load(file = "data/processed/cond_srat_beforeFiltering.RData")

UMAP with all cell types

plot_df = allcells_css@meta.data
plot_df = cbind(plot_df, allcells_css@reductions$umap_css@cell.embeddings)
plot_df$Condition = factor(plot_df$Condition, levels = c("healthy", "embolised", "regenerating"))
plot_df$allcells_simp[grepl("interaction", plot_df$allcells_simp)] = "Doublets"
plot_df$allcells_simp[grepl("Hepatocytes ", plot_df$allcells_simp)] = "Hepatocytes"
plot_df$allcells_simp = factor(plot_df$allcells_simp, 
                               levels = c("Cholangiocytes", "LSEC periportal", "LSEC pericentral", 
                                          "Endothelial cells (non-LSEC)", "Hepatocytes", "Kupffer cells",
                                          "Macrophages", "cDCs", "pDCs", "ab-T cells", "gd-T cells",
                                          "B cells", "Plasmablasts",
                                          "Stellate cells", "Dividing cells", "Doublets"))

pltumap = ggplot(plot_df, aes(y = UMAPCSS_2, x = UMAPCSS_1, colour = allcells_simp))+
  geom_point(size = pointsize)+
  labs(colour = "Cell Type")+
  scale_colour_manual(values = colsallct)+
  guides(colour = guide_legend(override.aes = list(size = 2)))+
  theme(aspect.ratio = 1,
        legend.title = element_text(hjust = 0),
        axis.line = element_blank(),
        axis.text = element_blank())

pdf("figure_panels/fig2/umap_all_celltypes.pdf", height = 5, width = 6, useDingbats = F)
print(pltumap)
dev.off()
png("figure_panels/fig2/umap_all_celltypes.png", height = 325, width = 400, antialias = "subpixel")
print(pltumap)
dev.off()

Proportion of each cell type per condition

df_cnt = table(plot_df$allcells_simp, plot_df$Condition)[-c(15,16),]
df_cnt_hep = reshape2::melt(apply(df_cnt, 2, function(x) round(x/sum(x)*100, 1))[5,])
df_cnt_nohep = reshape2::melt(apply(df_cnt[-5,], 2, function(x) round(x/sum(x)*100, 1)))

df_cnt_all = rbind(data.frame("Var1" = "Hepatocytes", "Var2" = rownames(df_cnt_hep), 
                              "value" = df_cnt_hep$value),
                   df_cnt_nohep)

mat_cnt_all = reshape2::dcast(data = df_cnt_all, formula = Var1 ~ Var2, value.var = "value")
rownames(mat_cnt_all) = mat_cnt_all$Var1
mat_cnt_all = mat_cnt_all[,-1]
ctord = hclust(dist(mat_cnt_all[,c(2,1,3)]))$order
ctord = c(7, 4, 10, 9, 5, 14, 11, 8, 3, 12,1,6,2,13)

annot_df = unique(plot_df[,c("major_ct","allcells_simp")])
annot_df = data.frame(row.names = annot_df$allcells_simp, "ct" = annot_df$major_ct)
heatp = pheatmap::pheatmap(mat_cnt_all[ctord,c(2,1,3)], cluster_cols = F, cluster_rows = F, 
                           treeheight_row = F, display_numbers = T, gaps_row = c(1,1), 
                           number_color = c("black", "white")[as.integer(mat_cnt_all[ctord,c(2,1,3)]>50)+1],
                           color = colorRampPalette(brewer.pal(n = 9, name = "Blues"))(100),
                           fontsize_row = 8, fontsize_col = 8,
                           annotation_row = annot_df, annotation_colors = list(ct = colsmajor[-6]))

pdf("figure_panels/fig2/heatmap_all_celltypes_prop.pdf", height = 3.5, width = 5, useDingbats = F)
print(heatp)
dev.off()

UMAP with conditions

randomrows = sample(1:nrow(plot_df), size = nrow(plot_df), replace = F)
pltumapcond = ggplot(plot_df[randomrows,], aes(y = UMAPCSS_2, x = UMAPCSS_1, colour = Condition))+
  geom_point(size = pointsize)+
  labs(colour = "Condition")+
  scale_colour_manual(values = colcond)+
  guides(colour = guide_legend(override.aes = list(size = 2)))+
  theme(aspect.ratio = 1,
        legend.title = element_text(hjust = 0),
        axis.line = element_blank(),
        axis.text = element_blank())

pdf("figure_panels/fig2/umap_condition.pdf", height = 5, width = 6, useDingbats = F)
print(pltumapcond)
dev.off()
png("figure_panels/fig2/umap_condition.png", height = 325, width = 400, antialias = "subpixel")
print(pltumapcond)
dev.off()

Top GO terms and load data

# loading and preparing data
go_enr_list = readRDS(file = "results/cond_effect/go_enr_list.RDS")

go_enr_major = go_enr_list$major_ct
for(cc in names(go_enr_major)){ # adding the terms for subsampled hepatocytes
  for(ct in names(go_enr_major[[cc]])){
    if(ct=="Hepatocytes"){
      go_enr_major[[cc]][[ct]] = go_enr_list$major_ct_hep[[cc]]$Hepatocytes
    }
    go_enr_major[[cc]][[ct]] = go_enr_major[[cc]][[ct]][go_enr_major[[cc]][[ct]]$DB=="GO Term",]
  }
}


for(comp in names(go_enr_major)){
  for(ct in names(go_enr_major[[comp]])){
    subdir = "figure_panels/fig2/goterms_cond_full/"
    go_enr_major[[comp]][[ct]]$geneID = unlist(go_enr_major[[comp]][[ct]]$geneID)
    write.table(go_enr_major[[comp]][[ct]], 
                file = paste0(subdir, "GOTermBP_", comp, "_", ct, ".csv"), 
                col.names = T, row.names = F, quote = F, sep = ",")
  }
}

Plotting GO Terms mixed between conditions

terms_list = list()
for(ct in unique(names(go_enr_major$healthy_v_embolised))){
  emb_terms = go_enr_major$healthy_v_embolised[[ct]][,c(2,7,8,11)]
  emb_terms = emb_terms[emb_terms$cond!="healthy",]
  reg_terms = go_enr_major$healthy_v_regenerating[[ct]][,c(2,7,8,11)]
  reg_terms = reg_terms[reg_terms$cond!="healthy",]
  both_terms = merge(emb_terms, reg_terms, by = 1, all = T)[,-c(4,7)]
  colnames(both_terms)[2:5] = c("qval_embolized", "genes_embolized", 
                                "qval_regenerating", "genes_regenerating")
  both_terms$qval_embolized[is.na(both_terms$qval_embolized)] = 1
  both_terms$qval_regenerating[is.na(both_terms$qval_regenerating)] = 1
  both_terms$genes_embolized[is.na(both_terms$genes_embolized)] = ""
  both_terms$genes_regenerating[is.na(both_terms$genes_regenerating)] = ""
  
  both_terms$genes_all = unlist(lapply(strsplit(paste0(both_terms$genes_embolized, "/",
                                                       both_terms$genes_regenerating),"/"), 
                                       function(x) paste0(unique(x[x!=""]), collapse = "/")))
  
  both_terms = both_terms[both_terms$qval_embolized<=0.05 | both_terms$qval_regenerating<=0.05,]
  
  both_terms_top = getTopTermsPaired(both_terms, ncl = 12, nt = 1)
  terms_list[[ct]] = both_terms_top
}

plt_list = list()
for(ct in names(terms_list)){
  sum_df = terms_list[[ct]][,c(1:3,7)]
  sum_df$Description = breakStr(as.character(sum_df$Description), n = 30)
  sum_df$qval_max[sum_df$qval_embolized>sum_df$qval_regenerating] = 1-sum_df$qval_max[sum_df$qval_embolized>sum_df$qval_regenerating]
  plot_df = reshape2::melt(sum_df[,1:3])
  plot_df$Description = factor(plot_df$Description, levels = sum_df$Description[order(sum_df$qval_max)])
  plot_df$variable = unlist(lapply(strsplit(as.character(plot_df$variable), "_"), function(x) x[2]))
  
  plt_list[[ct]] = ggplot(plot_df, aes(x = -log10(value), y = Description, fill = variable))+
    geom_col(position = "dodge")+
    geom_vline(xintercept = -log10(0.05), linetype = "dashed")+
    scale_x_continuous(limits = c(0, max(-log10(plot_df$value)+0.3)), expand = c(0,0))+
    scale_fill_manual(values = colcond)+
    theme(legend.position = "none",
          axis.ticks = element_line())
}

for(n in names(plt_list)){
  pdf(paste0("figure_panels/fig2/goterms_cond/goterms_mixedGO_", n, ".pdf"), 
      height = 2.8, width = 3)
  print(plt_list[[n]])
  dev.off()
}

cowplot::plot_grid(plotlist = plt_list, ncol = 4, align = "hv")

Top enriched GO Terms - group similar terms and report the top 1 of each of 5 groups detected

goplt = list()
gotab = list()
for(comp in names(go_enr_major)){
  for(gr in names(go_enr_major[[comp]])){
    conds = strsplit(comp, "_v_")[[1]]
    for(cc in conds){
      sub_df = go_enr_major[[comp]][[gr]][go_enr_major[[comp]][[gr]]$DB=="GO Term" &
                                            go_enr_major[[comp]][[gr]]$qvalue<=0.05 &
                                            go_enr_major[[comp]][[gr]]$cond==cc,]
      if(nrow(sub_df)>0){
        sub_df = getTopTerms(sub_df, ncl = 5, nt = 1)
        sub_df = sub_df[order(sub_df$qvalue, decreasing = F),]
        l = ifelse(nrow(sub_df)>10, 10, nrow(sub_df))
        sub_df = sub_df[1:l,]
        sub_df$geneID = unlist(sub_df$geneID)
        sub_df$Description = breakStr(as.character(sub_df$Description), 30)
        sub_df$Description = factor(sub_df$Description, 
                                    levels = rev(as.character(sub_df$Description)))
        
        goplt[[paste0(comp, "__", gr, "__", cc)]] = ggplot(sub_df, 
                                                           aes(x = -log10(qvalue), y = Description))+
          geom_col()+
          scale_x_continuous(expand = c(0,0), limits = c(0, max(-log10(sub_df$qvalue))+0.25))+
          labs(y = "GO Term Description", title = paste0(cc, " vs ", conds[conds!=cc]))+
          theme(axis.ticks.x = element_line(),
                axis.text.y = element_text(hjust = 1, vjust = 0.5, lineheight = 0.75,
                                           size = 7),
                plot.title = element_text(size = 10, hjust = 0),
                plot.title.position = "plot")
        gotab[[paste0(comp, "__", gr, "__", cc)]] = sub_df
      }
    }
  }
}
for(n in names(goplt)){
  pdf(paste0("figure_panels/fig2/goterms_cond/goterms_filtGO_", n, ".pdf"), 
      height = 2, width = 3)
  print(goplt[[n]])
  dev.off()
  
  write.csv(gotab[[n]], file = paste0("figure_panels/fig2/goterms_cond/goterms_filtGO_", n, ".csv"),
            row.names = F, quote = F)
}

DE between conditions for cell types

filt_comps = readRDS(file = "results/cond_effect/all_filt_cond_comps.RDS")

for(comp in names(filt_comps$major_ct)){
  subdir = "figure_panels/fig2/DE_tables/"
  de_df = filt_comps$major_ct[[comp]]
  de_df = de_df[de_df$celltype!="Hepatocytes",]
  de_df = rbind(de_df, filt_comps$major_ct_hep[[comp]])
  de_df = de_df[de_df$p_val_adj<=0.05,]
  
  write.csv(de_df[order(de_df$avg_logFC, decreasing = T),], 
            file = paste0(subdir, "DEcond_major_ct_", comp, ".csv"), 
            col.names = T, row.names = F, quote = F)
}

for(g in names(filt_comps)){
  for(comp in names(filt_comps[[g]])){
    subdir = file.path("to_send", g)
    write.csv(filt_comps[[g]][[comp]], file = paste0(subdir, "/DEconditions_", comp, "_", g, ".csv"), 
              col.names = T, row.names = F, quote = F)
  }
}

Major cell types

genes_to_plot = list("Hepatocytes" = c("HAMP", "SAA1", "FGA", "G0S2", "TNFRSF1A",
                                       "CYP3A7", "MT1A", "DDX21", "IL1RAP", "MT1X",
                                       "MGLL", "MFSD2A", "APOA4", "IRF7", "CEBPA"), 
                      "Immune" = c("S100A9", "CD81", "FCGR1A", "GBP5", "RUNX3", 
                                   "CD300E", "MNDA", "IL18BP", "IL4I1", "WDFY4",
                                   "TRDC", "TRBC1", "IFNG", "CD160", "IFITM1"), 
                      "Cholangiocytes" = c("CCL2", "IRF1", "KRT19", "CCL20", "ICAM1",
                                           "MAP3K12", "EPSTI1", "P4HA1", "PEAK1", "SEC24A",
                                           "LCN2", "RAMP1", "ITGB4", "BIK", "CCL28"), 
                      "Mesenchymal" = c(), 
                      "Endothelial" = c("SOX18", "PLCG2", "TFPI2", "KLF2", "STC1",
                                        "CXXC5", "NPR3", "CDK6", "ADGRG6", "PDGFB",
                                        "RFPL1", "ITGB3", "ADAMTS4", "KLF13", "GPR182"))

plt_list_cond = list()
df_list_cond = list()
for(ct in unique(filt_comps$major_ct$healthy_v_embolised$celltype)){
  g = if(ct=="Hepatocytes") "major_ct_hep" else "major_ct" #use correct Hep DE
  
  # prepare data frame
  plt_reg = filt_comps[[g]]$healthy_v_regenerating[,c(1,3,4,5)]
  plt_reg = unique(plt_reg[plt_reg$celltype==ct,])
  plt_emb = filt_comps[[g]]$healthy_v_embolised[,c(1,3,4,5)]
  plt_emb = unique(plt_emb[plt_emb$celltype==ct,])
  rownames(plt_reg) = plt_reg$gene
  rownames(plt_emb) = plt_emb$gene
  
  # add FC/pval for genes that are not in common
  plt_df = merge(plt_emb, plt_reg, by = 0, all = T)
  rownames(plt_df) = plt_df[,1]
  plt_df = plt_df[,-c(1,3,6,7)]
  colnames(plt_df) = c("gene", "FC_emb", "pval_emb", "FC_reg", "pval_reg")
  plt_df$gene = rownames(plt_df)
  plt_df$FC_emb[is.na(plt_df$FC_emb)] = 0
  plt_df$FC_reg[is.na(plt_df$FC_reg)] = 0
  plt_df$pval_emb[is.na(plt_df$pval_emb)] = 1
  plt_df$pval_reg[is.na(plt_df$pval_reg)] = 1
  
  # condition labels
  plt_df$cond = ifelse(plt_df$FC_emb<=(-0.2) & plt_df$pval_emb<=0.05 & plt_df$FC_reg>=0, "embolised",
                       ifelse(plt_df$FC_reg<=(-0.2) & plt_df$pval_reg<=0.05 & plt_df$FC_emb>=0,
                              "regenerating",
                              ifelse(plt_df$FC_reg<=(-0.2) & plt_df$pval_reg<=0.05 & 
                                       plt_df$FC_emb<=(-0.2) & plt_df$pval_emb<=0.05, "both","other")))
  plt_df$cond = factor(plt_df$cond, levels = c("both", "embolised", "regenerating", "other"))
  plt_df = plt_df[order(plt_df$cond, decreasing = T),]
  
  # genes to plot
  b_g_plt = genes_to_plot[[ct]][1:5]
  emb_g_plt = genes_to_plot[[ct]][6:10]
  reg_g_plt = genes_to_plot[[ct]][11:15]
  
  cols = c("both" = "#C54635", "regenerating" = "salmon", "embolised" = "darkred", "other" = "grey85")
  plt = ggplot(plt_df, aes(x = FC_emb, y = FC_reg, colour = cond))+
    geom_vline(xintercept = 0)+
    geom_hline(yintercept = 0)+
    geom_point()+
    geom_label_repel(data = plt_df[plt_df$gene %in% c(emb_g_plt, reg_g_plt, b_g_plt),], 
                     mapping = aes(label = gene), show.legend = F, min.segment.length = 0, 
                     size = 2.75, label.padding = 0.15)+
    scale_colour_manual(values = cols)+
    labs(title = ct, x = "logFC(healthy vs embolised)", colour = "Condition",
         y = "logFC(healthy vs regenerating)")+
    theme_bw()+
    theme(axis.text = element_text(colour = "black"),
          aspect.ratio = 1)
  
  plt_list_cond[[ct]] = plt
}

pdf("figure_panels/fig2/DE_cond_major_celltypes.pdf", height = 5, width = 6, useDingbats = F)
print(plt_list_cond)
dev.off()
pdf("figure_panels/fig2/DE_cond_major_celltypes_noLeg.pdf", height = 5, width = 6, useDingbats = F)
plt_list_cond = lapply(plt_list_cond, function(x) x+theme(legend.position = "none"))
print(plt_list_cond)
dev.off()

All cell types

plt_list_cond = list()
df_list_cond = list()
for(ct in unique(filt_comps$cell_type$healthy_v_embolised$celltype)){
  g = if(ct=="Hepatocytes") "cell_type_hep" else "cell_type" #use correct Hep DE
  
  # prepare data frame
  plt_reg = filt_comps[[g]]$healthy_v_regenerating[,c(1,3,4,5)]
  plt_reg = unique(plt_reg[plt_reg$celltype==ct,])
  plt_emb = filt_comps[[g]]$healthy_v_embolised[,c(1,3,4,5)]
  plt_emb = unique(plt_emb[plt_emb$celltype==ct,])
  rownames(plt_reg) = plt_reg$gene
  rownames(plt_emb) = plt_emb$gene
  
  # add FC/pval for genes that are not in common
  plt_df = merge(plt_emb, plt_reg, by = 0, all = T)
  rownames(plt_df) = plt_df[,1]
  plt_df = plt_df[,-c(1,3,6,7)]
  colnames(plt_df) = c("gene", "FC_emb", "pval_emb", "FC_reg", "pval_reg")
  plt_df$gene = rownames(plt_df)
  plt_df$FC_emb[is.na(plt_df$FC_emb)] = 0
  plt_df$FC_reg[is.na(plt_df$FC_reg)] = 0
  plt_df$pval_emb[is.na(plt_df$pval_emb)] = 1
  plt_df$pval_reg[is.na(plt_df$pval_reg)] = 1
  
  # condition labels
  plt_df$cond = ifelse(plt_df$FC_emb<=(-0.2) & plt_df$pval_emb<=0.05 & plt_df$FC_reg>=0, "embolised",
                       ifelse(plt_df$FC_reg<=(-0.2) & plt_df$pval_reg<=0.05 & plt_df$FC_emb>=0,
                              "regenerating",
                              ifelse(plt_df$FC_reg<=(-0.2) & plt_df$pval_reg<=0.05 & 
                                       plt_df$FC_emb<=(-0.2) & plt_df$pval_emb<=0.05, "both","other")))
  plt_df$cond = factor(plt_df$cond, levels = c("both", "embolised", "regenerating", "other"))
  plt_df = plt_df[order(plt_df$cond, decreasing = T),]
  
  # genes to plot
  ord = order(plt_df$FC_emb+plt_df$FC_reg, decreasing = F)
  b_g_plt = plt_df$gene[ord][plt_df$cond[ord]=="both" & plt_df$pval_emb[ord]<=0.05 & plt_df$pval_reg[ord]<=0.05][1:5]
  ord = order(plt_df$FC_emb, decreasing = F)
  emb_g_plt = plt_df$gene[ord][plt_df$cond[ord]=="embolised" & plt_df$pval_emb[ord]<=0.05][1:5]
  ord = order(plt_df$FC_reg, decreasing = F)
  reg_g_plt = plt_df$gene[ord][plt_df$cond[ord]=="regenerating" & plt_df$pval_reg[ord]<=0.05][1:5]
  
  cols = c("both" = "#C54635", "regenerating" = "salmon", "embolised" = "darkred", "other" = "grey85")
  plt = ggplot(plt_df, aes(x = FC_emb, y = FC_reg, colour = cond))+
    geom_vline(xintercept = 0)+
    geom_hline(yintercept = 0)+
    geom_point()+
    geom_label_repel(data = plt_df[plt_df$gene %in% c(emb_g_plt, reg_g_plt, b_g_plt),], 
                     mapping = aes(label = gene), show.legend = F, min.segment.length = 0, 
                     size = 2.75, label.padding = 0.15)+
    scale_colour_manual(values = cols)+
    labs(title = ct, x = "logFC(healthy vs embolised)", colour = "Condition",
         y = "logFC(healthy vs regenerating)")+
    theme_bw()+
    theme(axis.text = element_text(colour = "black"),
          aspect.ratio = 1)
  
  plt_list_cond[[ct]] = plt
}

pdf("figure_panels/fig2/DE_cond_all_celltypes.pdf", height = 5, width = 6, useDingbats = F)
print(plt_list_cond)
dev.off()
pdf("figure_panels/fig2/DE_cond_all_celltypes_noLeg.pdf", height = 5, width = 6, useDingbats = F)
plt_list_cond = lapply(plt_list_cond, function(x) x+theme(legend.position = "none"))
print(plt_list_cond)
dev.off()

Count DE genes per cell type in conditions - major cell types

de_df = filt_comps$major_ct$healthy_v_embolised
de_df = de_df[de_df$p_val_adj<=0.05 & abs(de_df$avg_logFC)>=0.2,]
genes_up = tapply(de_df$gene[de_df$avg_logFC>0], de_df$celltype[de_df$avg_logFC>0], list)
genes_dn = tapply(de_df$gene[de_df$avg_logFC<0], de_df$celltype[de_df$avg_logFC<0], list)
de_hep = filt_comps$major_ct_hep$healthy_v_embolised
genes_up$Hepatocytes = de_hep$gene[de_hep$avg_logFC>=0.2 & de_hep$p_val_adj<=0.05]
genes_dn$Hepatocytes = de_hep$gene[de_hep$avg_logFC<=(-0.2) & de_hep$p_val_adj<=0.05]
de_df_g = rbind(reshape2::melt(genes_up), reshape2::melt(genes_dn))
de_df_g$dedir = c(rep("higher in healthy", length(unlist(genes_up))), 
                  rep("higher in embolised", length(unlist(genes_dn))))

de_df_r = filt_comps$major_ct$healthy_v_regenerating
de_df_r = de_df_r[de_df_r$p_val_adj<=0.05 & abs(de_df_r$avg_logFC)>=0.2,]
genes_up = tapply(de_df_r$gene[de_df_r$avg_logFC>0], de_df_r$celltype[de_df_r$avg_logFC>0], list)
genes_dn = tapply(de_df_r$gene[de_df_r$avg_logFC<0], de_df_r$celltype[de_df_r$avg_logFC<0], list)
de_hep = filt_comps$major_ct_hep$healthy_v_regenerating
genes_up$Hepatocytes = de_hep$gene[de_hep$avg_logFC>=0.2 & de_hep$p_val_adj<=0.05]
genes_dn$Hepatocytes = de_hep$gene[de_hep$avg_logFC<=(-0.2) & de_hep$p_val_adj<=0.05]
de_df_r_g = rbind(reshape2::melt(genes_up), reshape2::melt(genes_dn))
de_df_r_g$dedir = c(rep("higher in healthy", length(unlist(genes_up))), 
                    rep("higher in regenerating", length(unlist(genes_dn))))


counts_de = rbind(data.frame(table(de_df_g$L1, de_df_g$dedir)), data.frame(table(de_df_r_g$L1, de_df_r_g$dedir)))
counts_de$comp = rep(c("healthy vs embolised", "healthy vs regenerating"), each = 10)
counts_de$Var2 = factor(counts_de$Var2, 
                        levels = c("higher in healthy", "higher in embolised", "higher in regenerating"))

plt = ggplot(counts_de, aes(x = Var1, y = Freq, fill = Var2))+
  facet_wrap(~comp, nrow = 2)+
  scale_fill_manual(values = c("higher in healthy" = "orange", "higher in regenerating" = "salmon", 
                               "higher in embolised" = "darkred"))+
  scale_y_continuous(expand = c(0,0), limits = c(0, 4000))+
  labs(fill = "DE genes", x = "Cell types", y = "# DE genes")+
  geom_bar(stat = "identity")+
  theme_classic()+
  theme(legend.justification=c(0,1), legend.position=c(0,1),
        legend.box.background = element_rect(colour = "transparent", fill = "transparent"),
        legend.background = element_rect(colour = "transparent", fill = "transparent"),
        legend.key.size = unit(0.4, "cm"),
        axis.text = element_text(colour = "black"),
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))

pdf("figure_panels/fig2/barplot_decond_major.pdf", height = 5, width = 3, useDingbats = F)
print(plt)
dev.off()
png("figure_panels/fig2/barplot_decond_major.png", height = 500, width = 350, antialias = "subpixel")
print(plt)
dev.off()

Count DE genes per cell type in conditions - all cell types

de_df = filt_comps$cell_type_simp$healthy_v_embolised
de_df = de_df[de_df$p_val_adj<=0.05 & abs(de_df$avg_logFC)>=0.2,]
genes_up = tapply(de_df$gene[de_df$avg_logFC>0], de_df$celltype[de_df$avg_logFC>0], list)
genes_dn = tapply(de_df$gene[de_df$avg_logFC<0], de_df$celltype[de_df$avg_logFC<0], list)
de_hep = filt_comps$major_ct_hep$healthy_v_embolised
genes_up$Hepatocytes = de_hep$gene[de_hep$avg_logFC>=0.2 & de_hep$p_val_adj<=0.05]
genes_dn$Hepatocytes = de_hep$gene[de_hep$avg_logFC<=(-0.2) & de_hep$p_val_adj<=0.05]
de_df_g = rbind(reshape2::melt(genes_up), reshape2::melt(genes_dn))
de_df_g$dedir = c(rep("higher in healthy", length(unlist(genes_up))), 
                  rep("higher in embolised", length(unlist(genes_dn))))

de_df_r = filt_comps$cell_type_simp$healthy_v_regenerating
de_df_r = de_df_r[de_df_r$p_val_adj<=0.05 & abs(de_df_r$avg_logFC)>=0.2,]
genes_up = tapply(de_df_r$gene[de_df_r$avg_logFC>0], de_df_r$celltype[de_df_r$avg_logFC>0], list)
genes_dn = tapply(de_df_r$gene[de_df_r$avg_logFC<0], de_df_r$celltype[de_df_r$avg_logFC<0], list)
de_hep = filt_comps$major_ct_hep$healthy_v_regenerating
genes_up$Hepatocytes = de_hep$gene[de_hep$avg_logFC>=0.2 & de_hep$p_val_adj<=0.05]
genes_dn$Hepatocytes = de_hep$gene[de_hep$avg_logFC<=(-0.2) & de_hep$p_val_adj<=0.05]
de_df_r_g = rbind(reshape2::melt(genes_up), reshape2::melt(genes_dn))
de_df_r_g$dedir = c(rep("higher in healthy", length(unlist(genes_up))), 
                    rep("higher in regenerating", length(unlist(genes_dn))))

counts_de = rbind(data.frame(table(de_df_g$L1, de_df_g$dedir)), 
                  data.frame(table(de_df_r_g$L1, de_df_r_g$dedir)))
counts_de$comp = rep(c("healthy vs embolised", "healthy vs regenerating"), each = 14)
counts_de$Var2 = factor(counts_de$Var2, 
                        levels = c("higher in healthy", "higher in embolised", "higher in regenerating"))

plt = ggplot(counts_de, aes(x = Var1, y = Freq, fill = Var2))+
  facet_wrap(~comp)+
  scale_fill_manual(values = c("higher in healthy" = "orange", 
                               "higher in regenerating" = "salmon", "higher in embolised" = "darkred"))+
  scale_y_continuous(expand = c(0,0), limits = c(0, 3500))+
  labs(fill = "DE genes", x = "Cell types", y = "# DE genes")+
  geom_bar(stat = "identity")+
  theme_classic()+
  theme(legend.justification=c(0,1), legend.position=c(0,1),
        legend.box.background = element_rect(colour = "transparent", fill = "transparent"),
        legend.background = element_rect(colour = "transparent", fill = "transparent"),
        legend.key.size = unit(0.4, "cm"),
        axis.text = element_text(colour = "black"),
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))

pdf("figure_panels/fig2/barplot_decond_all.pdf", height = 4, width = 6, useDingbats = F)
print(plt)
dev.off()
png("figure_panels/fig2/barplot_decond_all.png", height = 350, width = 500, antialias = "subpixel")
print(plt)
dev.off()

Supplementary

Gene and UMI counts per major cell type and condition

plot_df = allcells_css@meta.data
plot_df$Condition = factor(plot_df$Condition, levels = c("healthy", "embolised", "regenerating"))
plot_df$major_ct = factor(plot_df$major_ct, levels = c("Hepatocytes", "Endothelial", "Cholangiocytes", "Immune", "Mesenchymal"))

vio_counts_ct = ggplot(plot_df[!is.na(plot_df$major_ct),], 
                       aes(y = nCount_SCT, x = Condition, fill = Condition))+
  facet_wrap(~major_ct, scales = "free_x")+
  geom_violin()+
  scale_y_log10(name = "# UMI (normalised)")+
  scale_fill_manual(values = colcond)+
  theme(panel.grid.major.y = element_line(colour = "grey85"),
        legend.position = "none",
        axis.text = element_text(),
        axis.text.x = element_text(size = 9),
        axis.title.x = element_text(size = 10),
        strip.text = element_text(size = 9),
        strip.background = element_rect(fill = "transparent", colour = "black", size = 0.9))

pdf("figure_panels/fig2/violin_ncounts_major_celltype.pdf", height = 4, width = 7, useDingbats = F)
print(vio_counts_ct)
dev.off()
png("figure_panels/fig2/violin_ncounts_major_celltype.png", height = 425, width = 550, antialias = "subpixel")
print(vio_counts_ct)
dev.off()

vio_counts_ct = ggplot(plot_df[!is.na(plot_df$major_ct),], 
                       aes(y = nFeature_SCT, x = Condition, fill = Condition))+
  facet_wrap(~major_ct, scales = "free_x")+
  geom_violin()+
  scale_y_log10(name = "# Genes")+
  scale_fill_manual(values = colcond)+
  theme(panel.grid.major.y = element_line(colour = "grey85"),
        legend.position = "none",
        axis.text = element_text(),
        axis.text.x = element_text(size = 9),
        axis.title.x = element_text(size = 10),
        strip.text = element_text(size = 9),
        strip.background = element_rect(fill = "transparent", colour = "black", size = 0.9))

pdf("figure_panels/fig2/violin_ngenes_major_celltype.pdf", height = 4, width = 7, useDingbats = F)
print(vio_counts_ct)
dev.off()
png("figure_panels/fig2/violin_ngenes_major_celltype.png", height = 425, width = 550, antialias = "subpixel")
print(vio_counts_ct)
dev.off()

nUMI in filtered and kept cells

c_meta = Reduce(rbind, lapply(cond_srat, function(x) x@meta.data))[,c(11,12,13,14,15,17,23)]
h_meta = Reduce(rbind, lapply(healthy_srat, function(x) x@meta.data))[,c(11,12,13,14,15,17,23)]
all_meta = rbind(h_meta, c_meta)

all_meta$keep_allfilt = factor(c("No", "Yes")[all_meta$keep_allfilt+1], levels = c("Yes", "No"))
all_meta$Name = factor(all_meta$Name, 
                       levels = paste0("sc_", c(paste0("H",1:3),paste0("E",c(1:3,5:7)),paste0("R",c(1:3,5:7)))))
levels(all_meta$Name) = paste0("sc_", c(paste0("H",1:3),paste0("E",1:6),paste0("R",1:6)))

n_cells = all_meta %>% count(Name, Fraction, keep_allfilt)

plt = ggplot(all_meta, aes(x = Fraction, y = nCount_SCT, fill = keep_allfilt, colour = keep_allfilt))+
  facet_wrap(~Name, scales = "free", nrow = 3)+
  geom_violin(scale = "width")+
  scale_colour_manual(values = c("Yes" = "grey20","No" = "grey80"))+
  scale_fill_manual(values = c("Yes" = "grey20","No" = "grey80"))+
  scale_y_log10()+
  labs(y = "nUMI", colour = "keep", fill = "keep")+
  theme(strip.background = element_rect(fill = "transparent", colour = "black", size = 1),
        strip.text = element_text(margin = margin(2,0,2,0), size = 9),
        axis.ticks = element_line(), 
        axis.text.y = element_text(size = 6.5))

pdf("figure_panels/fig2/filtering_nUMI.pdf", height = 3.3, width = 7.1, useDingbats = F)
print(plt)
dev.off()

UMAP with healthy donors split

plot_df = data.frame(allcells_css@reductions$umap_css@cell.embeddings)
plot_df$Donor = factor(allcells_css@meta.data$Name)
plot_df$Donor = plyr::revalue(plot_df$Donor, c("sc_E5" = "sc_E4", "sc_R5" = "sc_R4", 
                                               "sc_E6" = "sc_E5", "sc_R6" = "sc_R5",
                                               "sc_E7" = "sc_E6", "sc_R7" = "sc_R6"))
plot_df$Donor = factor(plot_df$Donor, levels = c("sc_H1", "sc_H2", "sc_H3", "sc_E1", "sc_R1", "sc_E2", 
                                                 "sc_R2", "sc_E3", "sc_R3", "sc_E4", "sc_R4", "sc_E5", 
                                                 "sc_R5", "sc_E6", "sc_R6"))
plot_df$Condition = factor(allcells_css@meta.data$Condition, 
                           levels = c("healthy", "embolised", "regenerating"))
plot_df$major_ct = allcells_css@meta.data$major_ct
plot_df$allcells_major = allcells_css@meta.data$allcells_major
plot_df$major_ct[is.na(plot_df$major_ct)] = plot_df$allcells_major[is.na(plot_df$major_ct)]
plot_df$major_ct[plot_df$major_ct=="Dividing cells"] = "Immune" # the detected dividing cells are mostly (if not all) immune 
plot_df$major_ct = factor(plot_df$major_ct, levels = c("Hepatocytes", "Endothelial","Immune",
                                                       "Cholangiocytes","Mesenchymal", "Doublets"))

plt = ggplot(plot_df[sample(1:nrow(plot_df), nrow(plot_df), replace = F),], 
             aes(x = UMAPCSS_1, y = UMAPCSS_2, colour = major_ct))+
  facet_wrap(~Donor, drop = T, nrow = 3, ncol = 5)+
  geom_point(size = 0.08)+
  guides(colour = guide_legend(override.aes = list(size = 3), title = "Donors"))+
  scale_colour_manual(values = colsmajor)+
  theme_classic()+
  th_gen+
  theme(axis.line = element_blank(),
        axis.ticks = element_blank(),
        axis.title = element_blank(),
        axis.text = element_blank(),
        strip.text.x = element_text(size = 8, margin = margin(0.06,0,0.06,0, "cm")),
        strip.background.x = element_rect(size = 0.5, colour = "transparent"),
        legend.title = element_text(hjust = 0),
        aspect.ratio = 1)

pdf("figure_panels/fig2/umapAll_fresh_donors_split.pdf", 
    useDingbats = F, height = 4.1, width = 6.4)
print(plt)
dev.off()
png("figure_panels/fig2/umapAll_fresh_donors_split.png", 
    height=8, width=15, unit="cm", res=600, antialias = "subpixel")
print(plt)
dev.off()
pdf("figure_panels/fig2/umapAll_fresh_donors_noLeg_split.pdf", 
    useDingbats = F, height = 4.1, width = 5.8)
print(plt+theme(legend.position = "none"))
dev.off()
png("figure_panels/fig2/umapAll_fresh_donors_noLeg_split.png", 
    height=8, width=11, unit="cm", res=600, antialias = "subpixel")
print(plt+theme(legend.position = "none"))
dev.off()
pdf("figure_panels/fig2/umapAll_fresh_donors_split_cond.pdf", 
    useDingbats = F, height = 4.1, width = 6.4)
print(plt+facet_wrap(Condition~Donor, drop = T, nrow = 3, ncol = 5))
dev.off()
png("figure_panels/fig2/umapAll_fresh_donors_split_cond.png", 
    height=8, width=15, unit="cm", res=600, antialias = "subpixel")
print(plt+facet_wrap(Condition~Donor, drop = T, nrow = 3, ncol = 5))
dev.off()

Figure 3

Main Figure

Zonation between conditions

# gene expression
hep_cells = readRDS(file = "results/zonation_cond/hep_cells_zonation_rank.RDS")

# fitted expression
hep_sig_fits = readRDS(file = "results/zonation_cond/hep_sig_fits.RDS")
hep_fits_qval = readRDS(file = "results/zonation_cond/hep_fits_qval.RDS")

# GO and clusters
hep_res_list = readRDS(file = "results/zonation_cond/hep_got_clustering_list.RDS")

g_list = c( # embolised clusters
  "COX7C", # oxygen transport chain (also in regen 2)
  "CYP39A1", # steroid metabolism
  "ARG1", # organic acid catabolic process (also in regen 2)
  "PROX1", # developing liver marker
  "COPA", # vesicle trafficking (also in regen 1)
  "CDH2", # cluster 4
  "CD151", # cluster 4
           # regenerating clusters
  "AGPAT2", # lipid metabolism
  "COX6C", # electron transport chain
  "PTGR1", # fatty acid metabolism
  "FOXA3", # liver development (cluster 1)
  "FOXO1", # control of metabolism and CC arrest
  "CYP2A6", # cluster 4
  "CYP2B6", # cluster 4
  "CYP4A11", # cluster 4
  "SDC4", # cluster 4
  "STAT1", # cluster 4
  # normal zonation
  "SAA1", "SAA2", "HAMP", "C3", "CYP2E1","CYP1A2","HULC","BCHE","CYP3A4",
  "CRP", "SDS", "HAL", "IGFBP1", "IGFBP2", "BAAT", "SLCO1B3")

hep_df = data.frame(vals = c(hep_cells$healthy$zonation_pt, hep_cells$embolised$zonation_pt,
                              hep_cells$regenerating$zonation_pt), 
                     Condition = c(rep("healthy", length(hep_cells$healthy$zonation_pt)), 
                                   rep("embolised", length(hep_cells$embolised$zonation_pt)), 
                                   rep("regenerating", length(hep_cells$regenerating$zonation_pt))),
                     Donor = c(hep_cells$healthy$Donor, 
                               hep_cells$embolised$Donor, hep_cells$regenerating$Donor),
                     ct = c(hep_cells$healthy$allcells_simp, 
                               hep_cells$embolised$allcells_simp, hep_cells$regenerating$allcells_simp))
hep_df$bins10 = cut(hep_df$vals, 10)
hep_df$Condition = factor(hep_df$Condition, levels = c("healthy", "embolised", "regenerating"))
levels(hep_df$Donor) = c("sc_E1/R1", "sc_E2/R2", "sc_E3/R3", "sc_E4/R4", "sc_E5/R5", "sc_E6/R6", 
                         "sc_H1", "sc_H2", "sc_H3")

plt_dists = ggplot(hep_df, aes(x = bins10, fill = Condition))+
     facet_wrap(~Condition)+
     geom_bar()+
     labs(x = "zonation (binned)", y = "Number of cells")+
     scale_y_continuous(expand = c(0,0))+
     scale_fill_manual(values = colcond)+
     coord_flip()+
     guides(fill = guide_legend(title.position = "left"))+
     theme(axis.text.y = element_blank(),
           legend.title.align = 0,
           legend.position = "bottom",
           legend.box.margin = margin(0,0,0,0),
           legend.box.spacing = unit(0.05, "cm"))

pdf("figure_panels/fig3/hep_zonation_conditionDist.pdf", height = 2.3, width = 3.6, useDingbats = F)
print(plt_dists)
dev.off()

Donor density per condition

hep_don = ggplot(hep_df, aes(x = vals, colour = Donor))+
     facet_wrap(~Condition, scales = "free_y")+
     geom_hline(yintercept = 0)+
     geom_density()+
     labs(x = "zonation", y = "cell density")+
     scale_y_continuous(expand = c(0,0))+
     coord_flip()+
     guides(colour = guide_legend(title.position = "left", nrow = 3))+
     theme(axis.text.y = element_blank(),
           axis.ticks.x = element_line(),
           strip.text = element_text(margin = margin(0.04,0,0.04,0, "cm"), size = 7.5),
           legend.title.align = 0,
           legend.position = "bottom",
           legend.box.margin = margin(0,0,0,0),
           legend.box.spacing = unit(0.05, "cm"))

pdf("figure_panels/fig3/hep_zonation_conditionDonors.pdf", height = 3, width = 3.6, useDingbats = F)
print(hep_don)
dev.off()

Zonation grouped by donors

hep_don = ggplot(hep_df, aes(x = vals, colour = Condition))+
     facet_wrap(~Donor, scales = "free_y")+
     geom_hline(yintercept = 0)+
     geom_density()+
     labs(x = "zonation", y = "cell density")+
     scale_colour_manual(values = colcond)+
     scale_y_continuous(expand = c(0,0))+
     coord_flip()+
     guides(colour = guide_legend(title.position = "left", ncol = 3))+
     theme(axis.text.y = element_blank(),
           axis.ticks.x = element_line(),
           strip.text = element_text(margin = margin(0.04,0,0.04,0, "cm"), size = 7.5),
           legend.title.align = 0,
           legend.position = "bottom",
           legend.box.margin = margin(0,0,0,0),
           legend.box.spacing = unit(0.05, "cm"))

pdf("figure_panels/fig3/hep_zonation_DonorsCond.pdf", height = 3, width = 3.6, useDingbats = F)
print(hep_don)
dev.off()

Hep binned zonation normalised to healthy

hep_df = data.frame(vals = c(hep_cells$healthy$zonation_pt, hep_cells$embolised$zonation_pt,
                              hep_cells$regenerating$zonation_pt), 
                     Condition = c(rep("healthy", length(hep_cells$healthy$zonation_pt)), 
                                   rep("embolised", length(hep_cells$embolised$zonation_pt)), 
                                   rep("regenerating", length(hep_cells$regenerating$zonation_pt))))
hep_df$bins10 = cut(hep_df$vals, 10)
hep_df$Condition = factor(hep_df$Condition, levels = c("healthy", "embolised", "regenerating"))
hep_df = data.frame(table(hep_df$Condition, hep_df$bins10))
hep_df = hep_df[order(hep_df$Var1),]
hep_df$Freqnorm = unlist(tapply(hep_df$Freq, hep_df$Var1, function(x) x/sum(x)))
hep_df$Freqnorm_h = hep_df$Freqnorm/mean(hep_df$Freqnorm[hep_df$Var1=="healthy"])

plt_dists = ggplot(hep_df, aes(x = Var2, fill = Var1, y = Freqnorm_h))+
     facet_wrap(~Var1)+
     geom_col()+
     labs(x = "zonation (binned)", y = "Proportion of cells (compared to healthy)")+
     scale_y_continuous(expand = c(0,0), labels = c("0", "0.5", "1", "1.5", "2", "2.5"))+
     scale_fill_manual(values = colcond)+
     coord_flip()+
     guides(fill = guide_legend(title.position = "left"))+
     theme(axis.text.y = element_blank(),
           axis.text.x = element_text(size = 6),
           legend.position = "none")

pdf("figure_panels/fig3/hep_zonation_conditionDist_norm.pdf", 
    height = 2, width = 3.8, useDingbats = F)
print(plt_dists)
dev.off()

Plot GO Terms (general)

go_theme = theme_bw()+
  theme(axis.text = element_text(colour = "black", size = 7.5),
        axis.title = element_text(size = 8),
        plot.title = element_text(size = 9),
        plot.title.position = "plot",
        plot.subtitle = element_text(size = 7.5),
        axis.text.y = element_text(vjust = 0.6, colour = c("grey45", "grey17")),
        panel.grid.major.y = element_blank())

gopltFunc = function(plot_df, tit = "", subtit = ""){
  plot_df$Description = breakStr(plot_df$Description, 30)
  plot_df$Description = factor(plot_df$Description, levels = rev(plot_df$Description))
  plot_df$fill = as.character(rep(1:2, ceiling(nrow(plot_df)/2)))[1:nrow(plot_df)]
  plt = ggplot(plot_df, aes(x = -log10(qvalue), y = Description, fill = fill))+
    geom_col(show.legend = F)+
    scale_x_continuous(expand = c(0,0), lim = c(0, max(-log10(plot_df$qvalue))+1))+
    scale_fill_manual(values = c("grey17", "grey45"))+
    labs(title = tit, subtitle = subtit)+
    go_theme
  return(plt)
}

pdf("figure_panels/fig3/hep_go_emb_notcorr.pdf", height = 3.2, width = 3.6, useDingbats = F)
print(gopltFunc(hep_res_list$embolised$go_cl$all, tit = "GO Terms - Embolised", subtit = "cor < 0.3"))
dev.off()

pdf("figure_panels/fig3/hep_go_emb_corr.pdf", height = 3.2, width = 3.6, useDingbats = F)
print(gopltFunc(hep_res_list$embolised$go_cl$pos, tit = "GO Terms - Embolised", subtit = "cor < 0.3"))
dev.off()

pdf("figure_panels/fig3/hep_go_reg_notcorr.pdf", height = 3.2, width = 3.6, useDingbats = F)
print(gopltFunc(hep_res_list$regenerating$go_cl$all, tit = "GO Terms - Regenerating", subtit = "cor >= 0.3"))
dev.off()

pdf("figure_panels/fig3/hep_go_reg_corr.pdf", height = 3.2, width = 3.6, useDingbats = F)
print(gopltFunc(hep_res_list$regenerating$go_cl$pos, tit = "GO Terms - Regenerating", subtit = "cor >= 0.3"))
dev.off()

Plot GO Terms (each cluster)

pathgo = "figure_panels/fig3/goterms_clusters/"
dir.create(pathgo, showWarnings = FALSE)

for(cc in names(hep_res_list)){
  tit = ifelse(cc=="embolised", "GO Terms - Embolised", "GO Terms - Regenerating")
  for(n in 1:5){
    plot_df = hep_res_list[[cc]]$go_cl[[as.character(n)]]
    if(nrow(plot_df)>0){
      plt = gopltFunc(plot_df, tit = tit, subtit = paste0("Cluster ", n, "; cor < 0.3"))
      if(nrow(plot_df)%%2==1){
        plt = plt + theme(axis.text.y = element_text(vjust = 0.6, colour = c("grey17","grey45")))
      }
      pdf(paste0(pathgo, "hep_go_", cc, "_cl", n, ".pdf"), height = 3.2, width = 3.6, useDingbats = F)
      print(plt)
      dev.off()
    }
  }
}

Plots for individual genes

geneplot_theme = theme_bw()+
    theme(aspect.ratio = 1/1.5,
          panel.grid = element_blank(),
          axis.text = element_text(colour = "black", size = 7),
          axis.title = element_text(size = 8),
          plot.title = element_text(size = 10, face = "italic"),
          plot.title.position = "plot",
          legend.position = "bottom",
          legend.box.margin = margin(0,0,0,0),
          legend.box.spacing = unit(0.01,"cm"),
          legend.text = element_text(size = 8),
          legend.title = element_text(size = 9))

plt_gene_list = list()
for(g in g_list){
  plt_df = data.frame("exp" = c(hep_sig_fits$healthy[,g], hep_sig_fits$embolised[,g]),
                    "cond" = rep(c("healthy", "embolised"), each = 100),
                    "t" = rep(1:100, 2))
  plt_df$cond = factor(plt_df$cond, levels = c("healthy", "embolised"))
  pltemb = ggplot(plt_df, aes(x = t, y = exp, group = cond, col = cond))+
    geom_line(size = 1.5)+
    labs(x = "zonation", y = "Expression", title = g)+
    scale_colour_manual(values = colcond)+
    geneplot_theme
  
  plt_df = data.frame("exp" = c(hep_sig_fits$healthy[,g], hep_sig_fits$regenerating[,g]),
                    "cond" = rep(c("healthy", "regenerating"), each = 100),
                    "t" = rep(1:100, 2))
  plt_df$cond = factor(plt_df$cond, levels = c("healthy", "regenerating"))
  pltreg = ggplot(plt_df, aes(x = t, y = exp, group = cond, col = cond))+
    geom_line(size = 1.5)+
    labs(x = "zonation", y = "Expression", title = g)+
    scale_colour_manual(values = colcond)+
    geneplot_theme
  print(cowplot::plot_grid(pltemb, pltreg, ncol = 2))
  
  plt_gene_list[[g]] = list("embolised" = pltemb, "regenerating" = pltreg)
}
for(g in names(plt_gene_list)){
  for(n in names(plt_gene_list[[g]])){
    pdf(paste0("figure_panels/fig3/hep_zon_", g, "_",  n, ".pdf"), 
        height = 2.2, width = 2.5, useDingbats = F)
    print(plt_gene_list[[g]][[n]])
    dev.off()
  }
}

Plot individual genes with all three conditions

plt_3_list = list()
for(g in g_list){
  plt_df = data.frame("exp" = c(hep_sig_fits$healthy[,g], hep_sig_fits$embolised[,g],
                                hep_sig_fits$regenerating[,g]),
                    "cond" = rep(c("healthy", "embolised", "regenerating"), each = 100),
                    "t" = rep(1:100, 3))
  plt_df$cond = factor(plt_df$cond, levels = c("healthy", "embolised", "regenerating"))
  plt_3_list[[g]] = ggplot(plt_df, aes(x = t, y = exp, group = cond, col = cond))+
    geom_line(size = 1.5)+
    labs(x = "zonation", y = "Expression", title = g, colour = "Condition")+
    scale_colour_manual(values = colcond)+
    guides(colour = guide_legend(title.position = "top"))+
    geneplot_theme
}

for(g in names(plt_3_list)){
  pdf(paste0("figure_panels/fig3/hep_zon_", g, "_all.pdf"), 
      height = 2.5, width = 3, useDingbats = F)
  print(plt_3_list[[g]])
  dev.off()
}

Plot mean profiles

for(n in names(hep_res_list)){
  plot_df2 = hep_res_list[[n]]$df
  xxx = plot_df2[plot_df2$cond=="difference",]
  labsdf = data.frame(labs = c("higher in\nhealthy", paste0("higher in\n", n)), coord = c(0.5, -0.5))
  plt1 = ggplot(xxx, aes(x = Var1, y = value))+
      facet_grid(~cls, scales = "free")+
      geom_hline(yintercept = 0, colour = "grey60", size = 0.75, linetype = "dashed")+
      stat_summary(fun.data=mean_cl_boot, colour="grey40", alpha = 0.35, geom="linerange", group=1)+
      stat_summary(fun=mean, colour="black", geom="line", group=1)+
      labs(x = "zonation", y = "difference in\nscaled expression")+
      scale_y_continuous(limits = c(-1,1))+
      theme_bw()+
      theme(axis.text = element_text(colour = "black", size = 7),
            axis.title = element_text(size = 8))
  plt2 = ggplot()+
    geom_text(data = labsdf, mapping = aes(y = coord, label = labs), 
              x = 0.5, angle = 270, size = 3)+
    scale_x_continuous(limits = c(0,1))+
    scale_y_continuous(limits = c(-1,1))+
    cowplot::theme_nothing()
  plt = cowplot::plot_grid(plt1, plt2, ncol = 2, nrow = 1, rel_widths = c(1, 0.085), align = "hv")
  
  pdf(paste0("figure_panels/fig3/hep_zon_", n, "_mean.pdf"), height = 1.8, width = 6, useDingbats = F)
  print(plt)
  dev.off()
}

Plot heatmaps

# get profiles for plotting
usegenes = unique(c(as.character(hep_res_list$embolised$df$Var2),
                    as.character(hep_res_list$regenerating$df$Var2)))

bins_list = list()
for(n in names(hep_cells)){
  plot_df = cbind(data.frame("pt" = rep(1:10, each = 10)),hep_sig_fits[[n]][,usegenes])
  plot_df$bins10 = cut(plot_df$pt, 10)
  bins_list[[n]] = sapply(usegenes, function(x) scale(tapply(plot_df[,x], 
                                                             plot_df$bins10, mean))[10:1])
}
bins_mean = Reduce(rbind, bins_list)
colnames(bins_mean) = usegenes

# choose genes for plotting
ecl_genes = merge(unique(hep_res_list$embolised$df[,c(1,4)]), hep_fits_qval$embolised$fdr, 
                  by.x = 1,  by.y = 0)
top_e = ecl_genes %>% group_by(cls) %>% slice_min(order_by = pval_list, n = 5)
emb_heat = pheatmap(bins_mean[1:20,as.character(top_e$Var2)], show_rownames = F, 
                    fontsize_col = 6.2, border_color = NA, gaps_row = c(10,20), gaps_col = seq(5, 20, 5), 
                    clustering_method = "ward.D2", cluster_rows = F, cluster_cols = F, fontsize = 6)
pdf("figure_panels/fig3/heatmap_zonation_emb.pdf", height=2.5, width=3.8, useDingbats = F)
print(emb_heat)
dev.off()
write.csv(ecl_genes, file = "figure_panels/fig3/hep_embolised_clusters_fdr.csv", row.names = F, quote = F)

rcl_genes = merge(unique(hep_res_list$regenerating$df[,c(1,4)]), hep_fits_qval$regenerating$fdr, 
                  by.x = 1,  by.y = 0)
top_r = rcl_genes %>% group_by(cls) %>% slice_min(order_by = pval_list, n = 5)
reg_heat = pheatmap(bins_mean[c(1:10, 21:30),as.character(top_r$Var2)], show_rownames = F, 
                    fontsize_col = 7.5, border_color = NA, gaps_row = c(10,20), gaps_col = seq(5, 20, 5), 
                    clustering_method = "ward.D2", cluster_rows = F, cluster_cols = F, fontsize = 6.7)
pdf("figure_panels/fig3/heatmap_zonation_reg.pdf", height=2.5, width=3.8, useDingbats = F)
print(reg_heat)
dev.off()
write.csv(rcl_genes, file = "figure_panels/fig3/hep_regen_clusters_fdr.csv", row.names = F, quote = F)

Figure 4

Main Figure

Load data

only_end_cells = readRDS(file = "results/endothelial/only_end_cells_zon.RDS")
all_cor_end = read.csv("results/endothelial/correlations_zonation_endothelial.csv", 
                       header = T, row.names = 1)
end_fits_qval = readRDS(file = "results/endothelial/end_fits_qval.RDS")
mksimp = read.csv("results/endothelial/markers_endo_subpop_simp.csv", header = T, row.names = 1)
mksimp = mksimp[order(mksimp$avg_logFC, decreasing = T),]

top_genes = readRDS(file = "results/cond_effect/top_genes_mk.RDS")

go_end = readRDS("results/endothelial/GOTerms_correlations.RDS")

endo_col = c("Periportal LSEC" = "aquamarine2",  "Midzonal LSEC" = "aquamarine3",
             "Pericentral LSEC" = "aquamarine4", "LSEC (fenestr.)" = "darkslategray1",
             "LSEC (remodelling)" = "cyan4", "LSEC (interferon)" = "darkslateblue",
             "LSEC (high MT 1)" = "blueviolet", "LSEC (high MT 2)" = "darkviolet",
             "LSEC (stress)" = "lightblue", "EC non-LSEC" = "forestgreen", 
             "Lymphatic EC" = "chartreuse3", "Cycling cells" = "grey20")

UMAPs

only_end_cells = readRDS(file = "results/endothelial/only_end_cells_zon.RDS")
all_cor_end = read.csv("results/endothelial/correlations_zonation_endothelial.csv", 
                       header = T, row.names = 1)
end_fits_qval = readRDS(file = "results/endothelial/end_fits_qval.RDS")
mksimp = read.csv("results/endothelial/markers_endo_subpop_simp.csv", header = T, row.names = 1)
mksimp = mksimp[order(mksimp$avg_logFC, decreasing = T),]

top_genes = readRDS(file = "results/cond_effect/top_genes_mk.RDS")

go_end = readRDS("results/endothelial/GOTerms_correlations.RDS")

Frequency tables

df_cnt = table(only_end_cells$endo_simp, only_end_cells$Condition)
df_cnt_perCond = reshape2::melt(apply(df_cnt, 2, function(x) round(x/sum(x)*100, 1)))
df_cnt_perCl = reshape2::melt(apply(df_cnt, 1, function(x) round(x/sum(x)*100, 1)))

mat_cnt_all = reshape2::dcast(data = df_cnt_perCl, formula = Var1 ~ Var2, value.var = "value")
rownames(mat_cnt_all) = mat_cnt_all$Var1
mat_cnt_all = t(mat_cnt_all[,-1])
ctord = hclust(dist(mat_cnt_all[,c(2,1,3)]))$order

heatp = pheatmap::pheatmap(mat_cnt_all[ctord,c(2,1,3)], cluster_cols = F, cluster_rows = F, 
                           treeheight_row = F, display_numbers = T, fontsize_number = 7, 
                           number_color = c("black", "white")[as.integer(mat_cnt_all[ctord,c(2,1,3)]>50)+1],
                           color = colorRampPalette(brewer.pal(n = 9, name = "Blues"))(100),
                           fontsize_row = 7.5, fontsize_col = 7.5, angle_col = 0)

pdf("figure_panels/fig_endo/cluster_proportion_cond.pdf", height=3.2, width=3.3, useDingbats = F)
print(heatp)
dev.off()

Markers

mkplot = c("MGP", "AQP1", "CLEC14A", "EDNRB", "CLEC1B", "CLEC4G", "PLVAP", 
           "RBP7", "ANGPT2", "CTGF", "ISG15", "IFIT3", "MRO", "EGR1","PTGDS", "INMT", 
           "PROX1", "CCL21", "TOP2A")
mkplot = c("MGP", "AQP1", "CLEC14A", "IGFBP7", "CLEC1B", "CLEC4G", "PLVAP", 
           "RBP7", "ANGPT2", "CTGF", "ISG15", "IFIT3", "MRO", "EGR1","PTGDS", "INMT", 
           "PROX1", "CCL21", "TOP2A")

expdf = reshape2::melt(cbind(data.frame(Matrix::t(only_end_cells@assays$SCT@data[mkplot,]),
                                        "endo_simp" = only_end_cells$endo_simp)))
expdf$endo_simp = factor(expdf$endo_simp, levels = names(endo_col))

vioplt = ggplot(expdf, aes(x = endo_simp, y = value, fill = endo_simp))+
  facet_grid(variable~., scales = "free_y")+
  geom_violin(scale = "width", colour = "grey82", size = 0.35)+
  scale_fill_manual(values = endo_col)+
  scale_y_continuous(breaks = c(0, 2, 4))+
  labs(y = "Expression")+
  theme(legend.position = "none",
        axis.text.x = element_text(angle = -30, hjust = 0, vjust = 0.9),
        axis.title.x = element_blank(),
        axis.text.y = element_text(size = 6.5),
        axis.ticks = element_line(),
        strip.text.y = element_text(angle = 0, size = 7, margin = margin(0,0,0,0)),
        strip.background = element_rect(colour = "white", fill = "white"))

pdf("figure_panels_rev1/fig_endo/violins_mk_2.pdf", height=5.8, width=3.4, useDingbats = F)
print(vioplt)
dev.off()

boxplt = ggplot(expdf, aes(x = endo_simp, y = value, fill = endo_simp))+
  facet_grid(variable~., scales = "free_y")+
  geom_boxplot(colour = "grey50", outlier.shape = NA, size = 0.25)+
  scale_fill_manual(values = endo_col)+
  scale_y_continuous(breaks = c(0, 2, 4))+
  labs(y = "Expression")+
  theme(legend.position = "none",
        axis.text.x = element_text(angle = -30, hjust = 0, vjust = 0.9),
        axis.title.x = element_blank(),
        axis.text.y = element_text(size = 6.5),
        axis.ticks = element_line(),
        strip.text.y = element_text(angle = 0, size = 7, margin = margin(0,0,0,0)),
        strip.background = element_rect(colour = "white", fill = "white"))

pdf("figure_panels_rev1/fig_endo/boxes_mk_2.pdf", height=5.8, width=3.4, useDingbats = F)
print(boxplt)
dev.off()

Binned zonation

mkplot = c("MGP", "AQP1", "CLEC14A", "IGFBP7", "CLEC1B", "CLEC4M", "PLVAP", 
           "RBP7", "ANGPT2", "CTGF", "ISG15", "IFIT3", "MRO", "EGR1","PTGDS", "INMT", 
           "PROX1", "CCL21", "TOP2A")

expdf = reshape2::melt(cbind(data.frame(Matrix::t(only_end_cells@assays$SCT@data[mkplot,]),
                                        "endo_simp" = only_end_cells$endo_simp)))

Using endo_simp as id variables

expdf$endo_simp = factor(expdf$endo_simp, levels = names(endo_col))

vioplt = ggplot(expdf, aes(x = endo_simp, y = value, fill = endo_simp))+
  facet_grid(variable~., scales = "free_y")+
  geom_violin(scale = "width", colour = "grey82", size = 0.35)+
  scale_fill_manual(values = endo_col)+
  scale_y_continuous(breaks = c(0, 2, 4))+
  labs(y = "Expression")+
  theme(legend.position = "none",
        axis.text.x = element_text(angle = -30, hjust = 0, vjust = 0.9),
        axis.title.x = element_blank(),
        axis.text.y = element_text(size = 6.5),
        axis.ticks = element_line(),
        strip.text.y = element_text(angle = 0, size = 7, margin = margin(0,0,0,0)),
        strip.background = element_rect(colour = "white", fill = "white"))

Donor distributions

end_df = only_end_cells@meta.data[only_end_cells$goodcl,]

end_df = data.frame(vals = end_df$zonation_pt, 
                    Condition = end_df$Condition,
                    Donor = end_df$Donor,
                    ct = end_df$endo_simp)
end_df$bins10 = cut(end_df$vals, 10)
end_df$Condition = factor(end_df$Condition, levels = c("healthy", "regenerating", "embolised"))
end_df$Donor = factor(end_df$Donor)
end_df$Donor = plyr::revalue(end_df$Donor, c("HD1" = "sc_H1", "HD2" = "sc_H2", "HD3" = "sc_H3", 
                                             "DD1" = "sc_R1/E1", "DD2" = "sc_R2/E2", "DD3" = "sc_R3/E3",
                                             "DD5" = "sc_R4/E4", "DD6" = "sc_R5/E5", "DD7" = "sc_R6/E6"))
end_df$Donor = factor(end_df$Donor, levels = levels(end_df$Donor)[c(7:9, 1:6)])

end_don = ggplot(end_df, aes(x = vals, colour = Condition))+
  facet_wrap(~Donor, scales = "free_y")+
  geom_hline(yintercept = 0)+
  geom_density(size = 1.05)+
  labs(x = "Zonation", y = "Cell density")+
  scale_colour_manual(values = colcond)+
  scale_y_continuous(expand = c(0,0))+
  coord_flip()+
  guides(colour = guide_legend(title.position = "left", ncol = 3))+
  theme(axis.text.y = element_blank(),
        axis.ticks.x = element_line(),
        strip.text = element_text(margin = margin(0.04,0,0.04,0, "cm"), size = 7.5),
        strip.background = element_rect(fill = "white", colour = "white"),
        legend.title.align = 0,
        aspect.ratio = 0.8,
        legend.position = "bottom",
        legend.box.margin = margin(0,0,0,0),
        legend.box.spacing = unit(0.05, "cm"))

pdf("figure_panels/fig_endo/end_zonation_DonorsCond.pdf", height = 3.8, width = 3.4, useDingbats = F)
print(end_don)
dev.off()

Predicted zonation for other populations

plot_df = only_end_cells@meta.data[only_end_cells@meta.data$endo_simp %in% names(endo_col)[1:9],]
plot_df$bins100 = cut(plot_df$zonation_pt, 10)
plot_df$Condition = factor(plot_df$Condition, levels = c("healthy", "embolised", "regenerating"))
plot_df$endo_simp = factor(plot_df$endo_simp, levels = names(endo_col))

plot_ptall = ggplot(plot_df, aes(y = zonation_pt, x = endo_simp, fill = Condition))+
  facet_wrap(~endo_simp, scales = "free_x", ncol = 3)+
  geom_violin()+
  labs(y = "zonation", x = "clusters")+
  scale_y_continuous(expand = c(0,0), breaks = c(0, 0.25, 0.5, 0.75, 1), limits = c(0,1))+
  scale_fill_manual(values = colcond)+
  guides(fill = guide_legend(title.position = "left"))+
  theme_bw()+
  theme(axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.text.y = element_text(colour = "black", size = 7),
        axis.title = element_text(size = 8),
        strip.background = element_rect(fill = "white", colour = "black"),
        strip.text = element_text(size = 8),
        legend.title = element_blank(),
        legend.key.size = unit(0.3, "cm"),
        legend.box.spacing = unit(0.01, "cm"),
        legend.position = "bottom")

pdf("figure_panels/fig_endo/zonation_allLSEC.pdf", height=3.5, width=3.85, useDingbats = F)
print(plot_ptall)
dev.off()

Heatmap

only_end_cells_split = SplitObject(only_end_cells[,only_end_cells$goodcl], split.by = "Condition")

# filter markers from other cell types (can contaminate)
mk_others = unique(unlist(top_genes$cell_type[!grepl("LSEC ", names(top_genes$cell_type))]))
cor_sub = all_cor_end[!(all_cor_end$gene %in% mk_others),]

topgenes = unique(c(cor_sub$gene[order(cor_sub$HvE, decreasing = T)][1:30],
                    cor_sub$gene[order(cor_sub$HvR, decreasing = T)][1:30]))
bottomgenes = unique(c(cor_sub$gene[order(cor_sub$HvE, decreasing = F)][1:35],
                    cor_sub$gene[order(cor_sub$HvR, decreasing = F)][1:35]))

usegenes = c(topgenes, bottomgenes)

bins_list = list()
for(n in names(only_end_cells_split)){
  plot_df = cbind(data.frame("pt" = rep(1:10, each = 10)),end_fits_qval[[n]]$fits[,usegenes])
  plot_df$bins10 = cut(plot_df$pt, 10)
  bins_list[[n]] = sapply(usegenes, function(x) scale(tapply(plot_df[,x], 
                                                             plot_df$bins10, mean))[10:1])
}
bins_mean = Reduce(rbind, bins_list)

hct = hclust(dist(t(bins_list$healthy[,topgenes])), method = "ward.D2")
hcb = hclust(dist(t(bins_mean[,bottomgenes])), method = "ward.D2")
ord = c(topgenes[hct$order], bottomgenes[hcb$order])
cutb = factor(cutree(hcb, 10)[hcb$order], levels = unique(cutree(hcb, 8)[hcb$order]))

gaps = cumsum(table(cutb))[-10]+length(topgenes)

heat_plt = pheatmap(bins_mean[,ord], show_rownames = F, fontsize_col = 6.2, border_color = NA,
                    gaps_row = c(10,20), gaps_col = c(length(topgenes), gaps),
                    clustering_method = "ward.D2", cluster_rows = F, cluster_cols = F)

pdf("figure_panels/fig_endo/heatmap_corZon.pdf", height=3, width=7.8, useDingbats = F)
print(heat_plt)
dev.off()

GO Terms

go_theme = theme_bw()+
  theme(axis.text = element_text(colour = "black", size = 7.5),
        axis.title = element_text(size = 8),
        plot.title = element_text(size = 9),
        plot.title.position = "plot",
        plot.subtitle = element_text(size = 7.5),
        axis.text.y = element_text(vjust = 0.6, colour = c("grey45", "grey17")),
        panel.grid.major.y = element_blank())

gopltFunc = function(plot_df, tit = "", subtit = ""){
  plot_df$Description = breakStr(plot_df$Description, 30)
  plot_df$Description = factor(plot_df$Description, levels = rev(plot_df$Description))
  plot_df$fill = as.character(rep(1:2, ceiling(nrow(plot_df)/2)))[1:nrow(plot_df)]
  plt = ggplot(plot_df, aes(x = -log10(qvalue), y = Description, fill = fill))+
    geom_col(show.legend = F)+
    scale_x_continuous(expand = c(0,0), lim = c(0, max(-log10(plot_df$qvalue))+1))+
    scale_fill_manual(values = c("grey17", "grey45"))+
    labs(title = tit, subtitle = subtit)+
    go_theme
  return(plt)
}

pdf("figure_panels/fig_endo/end_go_emb_corr.pdf", height = 3.2, width = 3.6, useDingbats = F)
print(gopltFunc(go_end$cor_he, tit = "GO Terms - Embolised", subtit = "cor >= 0.3"))
dev.off()

pdf("figure_panels/fig_endo/end_go_reg_corr.pdf", height = 3.2, width = 3.6, useDingbats = F)
print(gopltFunc(go_end$cor_hr, tit = "GO Terms - Regenerating", subtit = "cor >= 0.3"))
dev.off()

Figure 5

Main Figure

Total interactions per cell type in each condition

net_names_l = readRDS(file = "results/cell_comm/updt/count_net_list.RDS")
# list per condition with interactions, genes, mean
reform_list = readRDS(file = "results/cell_comm/updt/reform_list.RDS")
# list per cond with ligand or receptor expression for each interaction
scoring_mats = readRDS(file = "results/cell_comm/updt/scoring_mats.RDS")
# reform_list + DE info + uniqueness/specificity of interaction/ligand/receptor
inter_df = readRDS(file = "results/cell_comm/updt/cond_diff_interact_DE.RDS")
# interactions qith certain functions involving which cell types and in which conditions
ct_g_cond_ann = readRDS(file = "results/cell_comm/updt/ct_g_cond_ann.RDS")
# interaction with L and R, pair of cell types, mean exp, and function
ct_int_exp_l = readRDS(file = "results/cell_comm/updt/interact_celltype_exp_group_list.RDS")
# annotation for all interactions
inter_annot = read.csv("data/interaction_annotation.csv", header = T)
inter_annot$funct[grepl("imm", inter_annot$funct)] = "immune"
inter_annot$funct[grepl("inflam", inter_annot$funct)] = "immune"
# mut info data
her_allint = readRDS(file = "./results/cell_comm/updt/interactions_mutInfo_condComp.RDS")
# GSEA on mutual info results
all_gsea_list = readRDS(file = "./results/cell_comm/updt/allgsea_interactions_mutInfo_condComp.RDS")
simp_gsea_list = readRDS(file = "./results/cell_comm/updt/simpgsea_interactions_mutInfo_condComp.RDS")
# network data
load("results/cell_comm/updt/median_networks_cond.RData")
load("results/cell_comm/updt/networks_cond.RData")
load("results/cell_comm/updt/median_networks_cond_midct.RData")
load("results/cell_comm/updt/median_networks_cond_umap.RData")
load("results/cell_comm/updt/networks_cond_umap.RData")
load("results/cell_comm/updt/median_networks_cond_umap_midct.RData")
comph_inters = readRDS(file = "results/cell_comm/updt/comph_inters.RDS")

Plot changes in number of interactions per condition

net_names_l = readRDS(file = "results/cell_comm/updt/count_net_list.RDS")

# list per condition with interactions, genes, mean
reform_list = readRDS(file = "results/cell_comm/updt/reform_list.RDS")
# list per cond with ligand or receptor expression for each interaction
scoring_mats = readRDS(file = "results/cell_comm/updt/scoring_mats.RDS")
# reform_list + DE info + uniqueness/specificity of interaction/ligand/receptor
inter_df = readRDS(file = "results/cell_comm/updt/cond_diff_interact_DE.RDS")
# interactions qith certain functions involving which cell types and in which conditions
ct_g_cond_ann = readRDS(file = "results/cell_comm/updt/ct_g_cond_ann.RDS")
# interaction with L and R, pair of cell types, mean exp, and function
ct_int_exp_l = readRDS(file = "results/cell_comm/updt/interact_celltype_exp_group_list.RDS")
# annotation for all interactions
inter_annot = read.csv("data/interaction_annotation.csv", header = T)
inter_annot$funct[grepl("imm", inter_annot$funct)] = "immune"
inter_annot$funct[grepl("inflam", inter_annot$funct)] = "immune"
# mut info data
her_allint = readRDS(file = "./results/cell_comm/updt/interactions_mutInfo_condComp.RDS")
# GSEA on mutual info results
all_gsea_list = readRDS(file = "./results/cell_comm/updt/allgsea_interactions_mutInfo_condComp.RDS")
simp_gsea_list = readRDS(file = "./results/cell_comm/updt/simpgsea_interactions_mutInfo_condComp.RDS")
# network data
load("results/cell_comm/updt/median_networks_cond.RData")
load("results/cell_comm/updt/networks_cond.RData")
load("results/cell_comm/updt/median_networks_cond_midct.RData")
load("results/cell_comm/updt/median_networks_cond_umap.RData")
load("results/cell_comm/updt/networks_cond_umap.RData")
load("results/cell_comm/updt/median_networks_cond_umap_midct.RData")
comph_inters = readRDS(file = "results/cell_comm/updt/comph_inters.RDS")

Non-ubiquitous interaction functions per cell type and condition

# correct some names
reform_list_int = list()
for(n in names(reform_list)[1:3]){
  reform_list[[n]]$ct1[grepl("Hepa", reform_list[[n]]$ct1)] = "Hepatocytes"
  reform_list[[n]]$ct2[grepl("Hepa", reform_list[[n]]$ct2)] = "Hepatocytes"
  reform_list[[n]]$ct1[grepl("LSEC_", reform_list[[n]]$ct1)] = "LSEC"
  reform_list[[n]]$ct2[grepl("LSEC_", reform_list[[n]]$ct2)] = "LSEC"
  reform_list_int[[n]] = reform_list[[n]][reform_list[[n]]$ct1!="Dividing cells" &
                                            reform_list[[n]]$ct2!="Dividing cells",]
  reform_list_int[[n]] = reform_list_int[[n]][!duplicated(reform_list_int[[n]][,1:3]),1:3]
}

# count interactions
infer_df = data.frame(expand.grid(unique(reform_list_int$healthy$ct1), 
                                  unique(reform_list_int$healthy$ct1)))
nr = nrow(infer_df)
infer_df = rbind(infer_df, infer_df, infer_df)
infer_df$count = 0
infer_df$cond = rep(names(reform_list)[1:3], each = nr)
colnames(infer_df)[1:2] = c("SOURCE", "TARGET")
for(i in 1:nrow(infer_df)){
  tmp = reform_list_int[[infer_df[i,4]]]
  infer_df[i,3] = dim(tmp[(tmp$ct1==infer_df[i,1] & tmp$ct2==infer_df[i,2]) |
                            (tmp$ct2==infer_df[i,1] & tmp$ct1==infer_df[i,2]),])[1]
}
## make one to count over all
infer_all = data.frame(expand.grid(unique(reform_list_int$healthy$ct1), 
                                   unique(reform_list_int$healthy$ct1)))
infer_all$count = 0
colnames(infer_all)[1:2] = c("SOURCE", "TARGET")
reform_all_int = unique(Reduce(rbind, reform_list_int))
for(i in 1:nrow(infer_all)){
  infer_all[i,3] = infer_all[i,3]+dim(tmp[(reform_all_int$ct1==infer_all[i,1] & reform_all_int$ct2==infer_all[i,2]) |
                                            (reform_all_int$ct2==infer_all[i,1] & reform_all_int$ct1==infer_all[i,2]),])[1]
}

infer_df_list = list(healthy = infer_df[infer_df$cond=="healthy",],
                     embolised = infer_df[infer_df$cond=="embolised",],
                     regenerating = infer_df[infer_df$cond=="regenerating",],
                     all = infer_all)

net_names_all = t(Reduce(rbind, lapply(infer_df_list[1:3], function(x) tapply(x$count, x$SOURCE, sum))))
colnames(net_names_all) = names(net_names_l)[1:3]
net_names_all = reshape2::melt(net_names_all)
net_names_all$Var1 = factor(net_names_all$Var1, 
                            levels = net_names_all[net_names_all$Var2=="healthy","Var1"][order(net_names_all[net_names_all$Var2=="healthy","value"], decreasing = F)])

total_unique = data.frame(t(t(Reduce(rbind, lapply(infer_df_list[4], 
                                                   function(x) tapply(x$count, x$SOURCE, sum))))))
total_unique$Var1 = rownames(total_unique)
colnames(total_unique)[1] = "value"
total_unique$Var2 = "total unique"

net_names_all$Var1 = factor(net_names_all$Var1, 
                            levels = total_unique$Var1[order(total_unique$value, 
                                                             decreasing = F)])
plt_counts = ggplot(net_names_all, aes(x = Var1, y = value))+
  geom_point(mapping = aes(colour = Var2), size = 1.6)+
  geom_point(data = total_unique, mapping = aes(shape = Var2), size = 1.6)+
  coord_flip()+
  scale_colour_manual(values = colcond)+
  scale_shape_manual(values = c(4))+
  scale_y_continuous(expand = c(0,0), limits = c(0,6500))+
  labs(y = "Total interactions", x = "Cell Type", colour = "Condition")+
  guides(colour = guide_legend(order = 1), shape = guide_legend(order = 2, title = NULL))+
  theme_bw()+
  theme(axis.text.x = element_text(angle = 30, vjust = 1, hjust = 1),
        axis.text = element_text(colour = "black"),
        axis.title.x = element_text(size = 8),
        legend.margin = margin(1, 1, 1, 1),
        legend.spacing.x = unit(0.05, "cm"),
        legend.direction = "horizontal",legend.box = "horizontal")
leg = cowplot::get_legend(plt_counts)

perc_df = data.frame("Var1" = net_names_all$Var1,
                     "Var2" = net_names_all$Var2,
                     "perc" = c(rep(0, length(unique(net_names_all$Var1))),
                                (net_names_all$value[net_names_all$Var2=="embolised"]-net_names_all$value[net_names_all$Var2=="healthy"])/net_names_all$value[net_names_all$Var2=="healthy"],
                                (net_names_all$value[net_names_all$Var2=="regenerating"]-net_names_all$value[net_names_all$Var2=="healthy"])/net_names_all$value[net_names_all$Var2=="healthy"])*100)
perc_df$col = ifelse(perc_df$perc>0, "green", "red")
perc_df = perc_df[perc_df$Var2!="healthy",]

perc_plt_emb = ggplot(perc_df[perc_df$Var2=="embolised",], aes(x = Var1, y = perc, fill = col))+
  geom_col()+
  scale_fill_manual(values = c("deepskyblue1", "red1"))+
  scale_y_continuous(expand = c(0,0), limits = c(-45,70))+
  coord_flip()+
  labs(y = "% change\nembolised vs healthy", x = "Cell Type")+
  theme_bw()+
  theme(axis.text.x = element_text(angle = 30, vjust = 1, hjust = 1, colour = "black"),
        axis.text.y = element_blank(),
        axis.title.x = element_text(size = 8),
        axis.title.y = element_blank(),
        axis.ticks.y = element_blank(),
        legend.position = "none")
perc_plt_reg = ggplot(perc_df[perc_df$Var2=="regenerating",], aes(x = Var1, y = perc, fill = col))+
  geom_col()+
  scale_fill_manual(values = c("deepskyblue1", "red1"))+
  scale_y_continuous(expand = c(0,0), limits = c(-45,70))+
  coord_flip()+
  labs(y = "% change\nregenerating vs healthy", x = "Cell Type")+
  theme_bw()+
  theme(axis.text.x = element_text(angle = 30, vjust = 1, hjust = 1, colour = "black"),
        axis.text.y = element_blank(),
        axis.title.x = element_text(size = 8),
        axis.title.y = element_blank(),
        axis.ticks.y = element_blank(),
        legend.position = "none")

cowplot::plot_grid(cowplot::plot_grid(plt_counts+theme(legend.position = "none"), 
                                      perc_plt_emb, perc_plt_reg, 
                                      ncol = 3, align = "h", axis = "r", rel_widths = c(1,0.5,0.5)),
                   leg, nrow = 2, rel_heights = c(1, 0.1))

pdf("figure_panels_rev1/fig_comm/mainFig_countsCond.pdf", height = 7, width = 7)

cowplot::plot_grid(cowplot::plot_grid(plt_counts+theme(legend.position = "none"), 
                                      perc_plt_reg, perc_plt_emb, 
                                      ncol = 3, align = "h", axis = "r", rel_widths = c(1,0.4,0.4)),
                   leg, nrow = 2, rel_heights = c(1, 0.1))
dev.off()

png 
  2 

Interaction functions per cell type and condition (using all interactions)

# interactions have to be absent in at least one condition
plot_df = ct_g_cond_ann[ct_g_cond_ann$celltypes %in% c("Stellate cells", "LSEC", "Kupffer cells", "cDCs",
                                                       "pDCs", "Endothelial cells (non-LSEC)") & 
                          ct_g_cond_ann$description %in% c("ECM", "inflammation", "migration",
                                                           "development", "immune activity", 
                                                           "immune suppression", "immune regulation"),]
plot_df$description = gsub(" ", "\n", plot_df$description, fixed = T)
plot_df$celltypes = gsub(" (", "\n(", plot_df$celltypes, fixed = T)
plot_df$description[grepl("immune", plot_df$description) | 
                      grepl("inflamm", plot_df$description)] = "immune"
bar_func_sub = ggplot(plot_df, aes(x = condition, fill = condition))+
  facet_grid(celltypes~description, scales = "free_y")+
  geom_bar()+
  scale_fill_manual(values = colcond)+
  theme_classic()+
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 7, colour = "black"),
        axis.text.y = element_text(size = 6.3, colour = "black"),
        axis.title = element_text(size = 7),
        axis.line = element_blank(),
        panel.background = element_rect(colour = "black"),
        panel.grid.major.y = element_line(),
        strip.text = element_text(size = 7),
        strip.background = element_blank(),
        aspect.ratio = 1,
        legend.position = "none")

pdf("figure_panels_rev1/fig_comm/mainFig_funcPart.pdf", height = 5.5, width = 4.5)
print(bar_func_sub)
dev.off()

null device 
          1 

plot_df = ct_g_cond_ann
plot_df$description = gsub(" ", "\n", plot_df$description, fixed = T)
plot_df$celltypes = gsub(" (", "\n(", plot_df$celltypes, fixed = T)
plot_df$description[grepl("immune", plot_df$description) | 
                      grepl("inflamm", plot_df$description)] = "immune"
bar_func_all = ggplot(plot_df, aes(x = condition, fill = condition))+
  facet_grid(celltypes~description, scales = "free_y")+
  geom_bar()+
  scale_fill_manual(values = colcond)+
  theme_classic()+
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 7, colour = "black"),
        axis.text.y = element_text(size = 6.3, colour = "black"),
        axis.title = element_text(size = 7),
        axis.line = element_blank(),
        panel.background = element_rect(colour = "black"),
        panel.grid.major.y = element_line(),
        strip.text = element_text(size = 7),
        strip.background = element_blank(),
        aspect.ratio = 1/2,
        legend.position = "none")

pdf("figure_panels_rev1/fig_comm/suppFig_countsFunc.pdf", paper = "a4", width = 9, height = 13.5)

print(bar_func_all)

*** recursive gc invocation
*** recursive gc invocation
*** recursive gc invocation
*** recursive gc invocation
*** recursive gc invocation
*** recursive gc invocation
*** recursive gc invocation
*** recursive gc invocation

dev.off()

null device 
          1 

GSEA qvalue for changing interactions (using all)

plot_df = Reduce(rbind, reform_list[1:3])
plot_df$cond = c(rep("healthy", nrow(reform_list$healthy)), 
                 rep("embolised", nrow(reform_list$embolised)), 
                 rep("regenerating", nrow(reform_list$regenerating)))
plot_df = unique(plot_df)
plot_df = data.frame(rbind(as.matrix(plot_df[,c(1,2,8)]), as.matrix(plot_df[,c(1,3,8)])))
plot_df = merge(plot_df, inter_annot, by = 1, all.x = T)
colnames(plot_df) = c("interaction", "celltypes", "condition", "description")

plot_df$description = gsub(" ", "\n", plot_df$description, fixed = T)
plot_df$celltypes = gsub(" (", "\n(", plot_df$celltypes, fixed = T)
plot_df$description[grepl("immune", plot_df$description) | 
                      grepl("inflamm", plot_df$description)] = "immune"
bar_func_all = ggplot(plot_df, aes(x = condition, fill = condition))+
  facet_grid(description~celltypes, scales = "free_y")+
  geom_bar()+
  scale_fill_manual(values = colcond)+
  theme_classic()+
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 7, colour = "black"),
        axis.text.y = element_text(size = 6.3, colour = "black"),
        axis.title = element_text(size = 7),
        axis.line = element_blank(),
        panel.background = element_rect(colour = "black"),
        panel.grid.major.y = element_line(),
        strip.text = element_text(size = 7),
        strip.background = element_blank(),
        aspect.ratio = 1,
        legend.position = "none")

Interaction expression matrices

plot_df = Reduce(rbind, all_gsea_list)
plot_df$q.val = -log10(plot_df$q.val+0.0000005)*(plot_df$sscore/abs(plot_df$sscore))
plot_df$comp = factor(plot_df$comp, levels = rev(c("he","hr","er")))
m = tapply(plot_df$q.val, plot_df$group,mean)
plot_df$group = factor(plot_df$group, levels = names(m)[order(m, decreasing = F)])
plot_df = plot_df[plot_df$comp!="er",]
plot_df$comp = as.character(plyr::revalue(plot_df$comp, 
                                          c("hr" = "regenerating", "he" = "embolised", "er" = "nothing")))
plot_df$comp = factor(plot_df$comp, levels = c("regenerating", "embolised"))

gseaplt = ggplot(plot_df, aes(x = q.val, y = group, fill = comp))+
  geom_col(position = "dodge")+
  geom_vline(xintercept = c(log10(0.05), -log10(0.05)), linetype = "dashed")+
  geom_vline(xintercept = c(0))+
  scale_fill_manual(values = colcond, drop = T)+
  guides(fill = guide_legend(title = "healthy vs", reverse = T))+
  labs(x = "q-value x score signal", y = "interaction type", fill = "comparison")+
  theme_bw()+
  theme(axis.text = element_text(colour = "black", size = 7),
        axis.title = element_text(colour = "black", size = 7.5),
        legend.position = c(0,1),
        legend.justification = c(-0.05,1.05),
        legend.title = element_text(size = 7.5),
        legend.text = element_text(size = 7),
        legend.margin = margin(0,0,0,0),
        legend.key.size = unit(0.35, "cm"))

pdf("figure_panels_rev1/fig_comm/mainFig_funcGSEA.pdf", height = 2.2, width = 2.85)
print(gseaplt)
dev.off()

null device 
          1 

Network plots - MDS

pltboth = ggplot()+
  geom_point(data = point_cond_df, mapping = aes(x = X1, y = X2, colour = ct2), 
             alpha = 0.6, show.legend = NA, size = 1.3)+
  geom_segment(data = pe_l[[2]], 
               mapping = aes(x = X1.x, xend = X1.y, y = X2.x, yend = X2.y, size = Freq), 
               alpha = 0.15, show.legend = F)+
  geom_point(data = pe_l[[1]], 
             mapping = aes(x = X1, y = X2, fill = ct), 
             alpha = 1, pch = 21, size = 4)+
  guides(colour = guide_legend(override.aes = list(size = 2.3)))+
  scale_size_continuous(range = c(0, 4), limits = range(pe_l[[2]]$Freq))+
  scale_colour_manual(values = colsmajor)+
  #scale_fill_manual(values = colsallct)+
  theme_classic()+ 
  theme(aspect.ratio = 1, 
        axis.line = element_blank(),
        axis.text = element_blank(),
        axis.title = element_blank(),
        axis.ticks = element_blank(),
        legend.title = element_blank(),
        legend.key.size = unit(0.35, "cm"),
        legend.spacing.y = unit(0.1, "cm"),
        legend.position = "none")
pdf("figure_panels_rev1/fig_comm/mainFig_inet_cond_all.pdf", 
    height = 3.1, width = 5.2, useDingbats = F)
print(pltboth)
dev.off()

plt_cond_l_ch = list()
for(cc in unique(pe_cond_l_ch[[2]]$condition)){
  plt_cond_l_ch[[cc]] = ggplot()+
    geom_segment(data = pe_cond_l_ch[[2]][pe_cond_l_ch[[2]]$condition==cc,], 
                 mapping = aes(x = X1.x, xend = X1.y, y = X2.x, yend = X2.y, size = Freq, alpha = Freq))+
    geom_point(data = pe_cond_l_ch[[1]], 
               mapping = aes(x = X1, y = X2, fill = ct), 
               alpha = 1, pch = 21, size = 4)+
    scale_size_continuous(range = c(0, 4), limits = range(pe_cond_l_ch[[2]]$Freq))+
    scale_alpha_continuous(range = c(0, 1), limits = range(pe_cond_l_ch[[2]]$Freq))+
    #scale_fill_manual(values = colsallct)+
    labs(title = cc)+
    guides(size = guide_legend(direction = "horizontal", nrow = 2),
           alpha = guide_legend(direction = "horizontal", nrow = 2))+
    theme_classic()+ 
    theme(aspect.ratio = 1, 
          axis.line = element_blank(),
          axis.text = element_blank(),
          axis.title = element_blank(),
          axis.ticks = element_blank(),
          panel.background = element_rect(colour = "black"),
          legend.title = element_blank(),
          legend.key.size = unit(0.35, "cm"),
          legend.position = "none")
}
plt_cond_l_ch[["leg"]] = cowplot::get_legend(plt_cond_l_ch[["healthy"]])

pdf("figure_panels_rev1/fig_comm/mainFig_inet_cond_median.pdf", height = 3.2, width = 5, useDingbats = F)
print(plt_cond_l_ch$healthy)
print(plt_cond_l_ch$embolised)
print(plt_cond_l_ch$regenerating)
dev.off()

Network plots with middle cell type resolution - MDS

pltboth = ggplot()+
  geom_point(data = point_cond_df, mapping = aes(x = X1, y = X2, colour = ct2), 
             alpha = 0.6, show.legend = NA, size = 1.3)+
  geom_segment(data = pe_l[[2]], 
               mapping = aes(x = X1.x, xend = X1.y, y = X2.x, yend = X2.y, size = Freq), 
               alpha = 0.15, show.legend = F)+
  geom_point(data = pe_l[[1]], 
             mapping = aes(x = X1, y = X2, fill = ct), 
             alpha = 1, pch = 21, size = 4)+
  guides(colour = guide_legend(override.aes = list(size = 2.3)))+
  scale_size_continuous(range = c(0, 4), limits = range(pe_l[[2]]$Freq))+
  scale_colour_manual(values = colsmajor)+
  #scale_fill_manual(values = colsallct)+
  theme_classic()+ 
  theme(aspect.ratio = 1, 
        axis.line = element_blank(),
        axis.text = element_blank(),
        axis.title = element_blank(),
        axis.ticks = element_blank(),
        legend.title = element_blank(),
        legend.key.size = unit(0.35, "cm"),
        legend.spacing.y = unit(0.1, "cm"),
        legend.position = "none")
pdf("figure_panels_rev1/fig_comm/mainFig_inet_cond_all.pdf", 
    height = 3.1, width = 5.2, useDingbats = F)
print(pltboth)
dev.off()

null device 
          1 

plt_cond_l_ch = list()
for(cc in unique(pe_cond_l_ch[[2]]$condition)){
  plt_cond_l_ch[[cc]] = ggplot()+
    geom_segment(data = pe_cond_l_ch[[2]][pe_cond_l_ch[[2]]$condition==cc,], 
                 mapping = aes(x = X1.x, xend = X1.y, y = X2.x, yend = X2.y, size = Freq, alpha = Freq))+
    geom_point(data = pe_cond_l_ch[[1]], 
               mapping = aes(x = X1, y = X2, fill = ct), 
               alpha = 1, pch = 21, size = 4)+
    scale_size_continuous(range = c(0, 4), limits = range(pe_cond_l_ch[[2]]$Freq))+
    scale_alpha_continuous(range = c(0, 1), limits = range(pe_cond_l_ch[[2]]$Freq))+
    #scale_fill_manual(values = colsallct)+
    labs(title = cc)+
    guides(size = guide_legend(direction = "horizontal", nrow = 2),
           alpha = guide_legend(direction = "horizontal", nrow = 2))+
    theme_classic()+ 
    theme(aspect.ratio = 1, 
          axis.line = element_blank(),
          axis.text = element_blank(),
          axis.title = element_blank(),
          axis.ticks = element_blank(),
          panel.background = element_rect(colour = "black"),
          legend.title = element_blank(),
          legend.key.size = unit(0.35, "cm"),
          legend.position = "none")
}
plt_cond_l_ch[["leg"]] = cowplot::get_legend(plt_cond_l_ch[["healthy"]])

pdf("figure_panels_rev1/fig_comm/mainFig_inet_cond_median.pdf", height = 3.2, width = 5, useDingbats = F)

print(plt_cond_l_ch$healthy)
print(plt_cond_l_ch$embolised)

print(plt_cond_l_ch$regenerating)
dev.off()

null device 
          1 

Network plots - UMAP

pltboth = ggplot()+
  geom_point(data = point_cond_umap_df, mapping = aes(x = X1, y = X2, colour = ct2), 
             alpha = 0.6, show.legend = NA, size = 1.3)+
  geom_segment(data = pe_umap_l[[2]], 
               mapping = aes(x = X1.x, xend = X1.y, y = X2.x, yend = X2.y, size = Freq), 
               alpha = 0.15, show.legend = F)+
  geom_point(data = pe_umap_l[[1]], 
             mapping = aes(x = X1, y = X2, fill = ct), 
             alpha = 1, pch = 21, size = 4)+
  guides(colour = guide_legend(override.aes = list(size = 2.3)))+
  scale_size_continuous(range = c(0, 4), limits = range(pe_l[[2]]$Freq))+
  scale_colour_manual(values = colsmajor)+
  scale_fill_manual(values = colsallct)+
  theme_classic()+ 
  theme(aspect.ratio = 1, 
        axis.line = element_blank(),
        axis.text = element_blank(),
        axis.title = element_blank(),
        axis.ticks = element_blank(),
        legend.title = element_blank(),
        legend.key.size = unit(0.35, "cm"),
        legend.spacing.y = unit(0.1, "cm"),
        legend.position = "none")
pdf("figure_panels_rev1/fig_comm/mainFig_inet_cond_umap_all.pdf", height = 3.1, width = 5.2, useDingbats = F)
print(pltboth)
dev.off()

plt_cond_l_ch = list()
for(cc in unique(pe_umap_cond_l_ch[[2]]$cond)){
  plt_cond_l_ch[[cc]] = ggplot()+
    geom_segment(data = pe_umap_cond_l_ch[[2]][pe_umap_cond_l_ch[[2]]$cond==cc,], 
                 mapping = aes(x = X1.x, xend = X1.y, y = X2.x, yend = X2.y, size = Freq, alpha = Freq))+
    geom_point(data = pe_umap_cond_l_ch[[1]], 
               mapping = aes(x = X1, y = X2, fill = ct), 
               alpha = 1, pch = 21, size = 4)+
    scale_size_continuous(range = c(0, 4), limits = range(pe_umap_cond_l_ch[[2]]$Freq))+
    scale_alpha_continuous(range = c(0, 1), limits = range(pe_umap_cond_l_ch[[2]]$Freq))+
    scale_fill_manual(values = colsallct)+
    labs(title = cc)+
    guides(size = guide_legend(direction = "horizontal", nrow = 2),
           alpha = guide_legend(direction = "horizontal", nrow = 2))+
    theme_classic()+ 
    theme(aspect.ratio = 1, 
          axis.line = element_blank(),
          axis.text = element_blank(),
          axis.title = element_blank(),
          axis.ticks = element_blank(),
          panel.background = element_rect(colour = "black"),
          legend.title = element_blank(),
          legend.key.size = unit(0.35, "cm"),
          legend.position = "none")
}
plt_cond_l_ch[["leg"]] = cowplot::get_legend(plt_cond_l_ch[["healthy"]])

pdf("figure_panels_rev1/fig_comm/mainFig_inet_cond_umap_median.pdf", height = 3.2, width = 5, useDingbats = F)
print(plt_cond_l_ch$healthy)
print(plt_cond_l_ch$embolised)
print(plt_cond_l_ch$regenerating)
dev.off()

Network plots with middle cell type resolution - UMAP

pltboth = ggplot()+
  geom_point(data = point_cond_umap_df, mapping = aes(x = X1, y = X2, colour = ct2), 
             alpha = 0.6, show.legend = NA, size = 1.3)+
  geom_segment(data = pe_umap_l[[2]], 
               mapping = aes(x = X1.x, xend = X1.y, y = X2.x, yend = X2.y, size = Freq), 
               alpha = 0.15, show.legend = F)+
  geom_point(data = pe_umap_l[[1]], 
             mapping = aes(x = X1, y = X2, fill = ct), 
             alpha = 1, pch = 21, size = 4)+
  guides(colour = guide_legend(override.aes = list(size = 2.3)))+
  scale_size_continuous(range = c(0, 4), limits = range(pe_l[[2]]$Freq))+
  scale_colour_manual(values = colsmajor)+
  scale_fill_manual(values = colsallct)+
  theme_classic()+ 
  theme(aspect.ratio = 1, 
        axis.line = element_blank(),
        axis.text = element_blank(),
        axis.title = element_blank(),
        axis.ticks = element_blank(),
        legend.title = element_blank(),
        legend.key.size = unit(0.35, "cm"),
        legend.spacing.y = unit(0.1, "cm"),
        legend.position = "none")
pdf("figure_panels_rev1/fig_comm/mainFig_inet_cond_umap_all.pdf", height = 3.1, width = 5.2, useDingbats = F)
print(pltboth)
dev.off()

null device 
          1 

plt_cond_l_ch = list()
for(cc in unique(pe_umap_cond_l_ch[[2]]$cond)){
  plt_cond_l_ch[[cc]] = ggplot()+
    geom_segment(data = pe_umap_cond_l_ch[[2]][pe_umap_cond_l_ch[[2]]$cond==cc,], 
                 mapping = aes(x = X1.x, xend = X1.y, y = X2.x, yend = X2.y, size = Freq, alpha = Freq))+
    geom_point(data = pe_umap_cond_l_ch[[1]], 
               mapping = aes(x = X1, y = X2, fill = ct), 
               alpha = 1, pch = 21, size = 4)+
    scale_size_continuous(range = c(0, 4), limits = range(pe_umap_cond_l_ch[[2]]$Freq))+
    scale_alpha_continuous(range = c(0, 1), limits = range(pe_umap_cond_l_ch[[2]]$Freq))+
    scale_fill_manual(values = colsallct)+
    labs(title = cc)+
    guides(size = guide_legend(direction = "horizontal", nrow = 2),
           alpha = guide_legend(direction = "horizontal", nrow = 2))+
    theme_classic()+ 
    theme(aspect.ratio = 1, 
          axis.line = element_blank(),
          axis.text = element_blank(),
          axis.title = element_blank(),
          axis.ticks = element_blank(),
          panel.background = element_rect(colour = "black"),
          legend.title = element_blank(),
          legend.key.size = unit(0.35, "cm"),
          legend.position = "none")
}
plt_cond_l_ch[["leg"]] = cowplot::get_legend(plt_cond_l_ch[["healthy"]])

pdf("figure_panels_rev1/fig_comm/mainFig_inet_cond_umap_median.pdf", height = 3.2, width = 5, useDingbats = F)

print(plt_cond_l_ch$healthy)
print(plt_cond_l_ch$embolised)

print(plt_cond_l_ch$regenerating)
dev.off()

null device 
          1 

Heatmaps for interactions with immune cells

# interactions unique to healthy or to emb/regen
comph_inters = c(setdiff(inter_df$healthy$id_cp_interaction,
                         c(inter_df$embolised$id_cp_interaction, inter_df$regenerating$id_cp_interaction)),
                 setdiff(inter_df$embolised$id_cp_interaction, inter_df$healthy$id_cp_interaction),
                 setdiff(inter_df$regenerating$id_cp_interaction, inter_df$healthy$id_cp_interaction))

edge_min = edge_cond_df[,c(1:4,9:12)]
edge_min = merge(edge_min, unique(her_allint_f[,c(1,13)]), by.x = 3, by.y = 1, all.x = T)
edge_min = merge(edge_min, unique(her_allint_f[,c(1,17)]), by.x = 1, by.y = 1, all.x = T)
edge_min = edge_min[edge_min$maj_g1=="Immune" | edge_min$maj_g2=="Immune",]
edge_min = edge_min[edge_min$ct_g1!=edge_min$ct_g2,]
edge_min = edge_min[edge_min$id_cp_interaction %in% comph_inters,]

devdf = filtFunc(her_allint_f[her_allint_f$funct=="development",], expthr = 0.15)
ecmdf = filtFunc(her_allint_f[her_allint_f$funct=="ECM",], expthr = 0.15)

plot_df = filtFunc(her_allint_f[her_allint_f$id_cp_interaction %in% unique(edge_min$id_cp_interaction),],
                   expthr = 0.07, minct = 2)
plot_df = plot_df[!(plot_df$intlr %in% c(devdf$intlr, ecmdf$intlr)),]
imm_plt = pltFunc(plot_df)
pdf("figure_panels_rev1/fig_comm/mainFig_immuneNetInteract_many.pdf", height = 3, width = 20)
print(imm_plt)
dev.off()

plot_df = filtFunc(her_allint_f[her_allint_f$id_cp_interaction %in% unique(edge_min$id_cp_interaction),],
                   expthr = 0.3, minct = 2)
plot_df = plot_df[!(plot_df$intlr %in% c(devdf$intlr, ecmdf$intlr)),]
imm_plt = pltFunc(plot_df)
pdf("figure_panels_rev1/fig_comm/mainFig_immuneNetInteract.pdf", height = 2.86, width = 7.1)
print(imm_plt)
dev.off()

Heatmaps with number of interactions per condition

cnts_list = list()
tab_list = list()
for(n in names(inter_df)[1:3]){
  subdfct = unique(inter_df[[n]][,1:3])
  subdfct$ct1[grepl("LSEC_", subdfct$ct1)] = "LSEC"
  subdfct$ct2[grepl("LSEC_", subdfct$ct2)] = "LSEC"
  subdfct$ct1[grepl("Hepat", subdfct$ct1)] = "Hepatocytes"
  subdfct$ct2[grepl("Hepat", subdfct$ct2)] = "Hepatocytes"
  #subdfct$ct1[grepl("non-", subdfct$ct1)] = "Endothelial cells\n(non-LSEC)"
  #subdfct$ct2[grepl("non-", subdfct$ct2)] = "Endothelial cells\n(non-LSEC)"
  subdfct = unique(subdfct)
  dfct = table(data.frame("ct1" = c(subdfct$ct1, subdfct$ct2),
                          "ct2" = c(subdfct$ct2, subdfct$ct1)))
  hcl = hclust(dist(dfct), method = "ward.D2")
  plot_df = data.frame(dfct)
  plot_df$ct1 = factor(plot_df$ct1, levels = levels(net_names_all$Var1))
  plot_df$ct2 = factor(plot_df$ct2, levels = levels(net_names_all$Var1))
  #plot_df$ct1 = factor(plot_df$ct1, levels = rev(c("Stellate cells", "LSEC", "Endothelial cells\n(non-LSEC)",
  #                                                 "Kupffer cells", "cDCs", "Macrophages", "Cholangiocytes",
  #                                                 "pDCs", "ab-T cells", "Hepatocytes", "gd-T cells", 
  #                                                 "B cells", "Plasmablasts")))
  #plot_df$ct2 = factor(plot_df$ct2, levels = rev(c("Stellate cells", "LSEC", "Endothelial cells\n(non-LSEC)",
  #                                                 "Kupffer cells", "cDCs", "Macrophages", "Cholangiocytes",
  #                                                 "pDCs", "ab-T cells", "Hepatocytes", "gd-T cells", 
  #                                                 "B cells", "Plasmablasts")))
  pltcnts = ggplot(plot_df, aes(x = ct1, y = ct2, size = Freq, colour = Freq))+
    geom_point()+
    labs(x = "cell type", y = "cell type", title = n)+
    guides(colour = guide_legend(title = "# interactions", reverse = T), 
           size = guide_legend(title = "# interactions", reverse = T))+
    scale_colour_gradientn(colors = c("#ededed", "#deebf7", "#9ecae1", "#bcbddc", "#756bb1"), 
                           limits = c(0, 240))+
    scale_size_continuous(limits = c(0, 240))+
    theme_classic()+
    theme(aspect.ratio = 1,
          panel.background = element_rect(colour = "black"),
          axis.line = element_blank(),
          axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.4),
          axis.text = element_text(colour = "black"))
  cnts_list[[n]] = pltcnts
  tab_list[[n]] = dfct
}
cnts_list$legend = cowplot::get_legend(cnts_list$healthy)
cnts_list$healthy = cnts_list$healthy+theme(legend.position = "none")
cnts_list$embolised = cnts_list$embolised+theme(legend.position = "none")
cnts_list$regenerating = cnts_list$regenerating+theme(legend.position = "none")

pdf("figure_panels_rev1/fig_comm/suppFig_countsCond.pdf", height = 15, width = 15)
cnts_list[[1]]+cnts_list[[3]]+cnts_list[[2]]+cnts_list[[4]]+      patchwork::plot_layout(ncol = 2)
dev.off()

Heatmaps with difference in number of interactions vs healthy

Heatmaps by mid cell type

diff_cnt_l = list()
for(n in c("embolised", "regenerating")){
  dif_df = tab_list[[n]]-tab_list$healthy
  hcl = hclust(dist(dif_df), method = "ward.D2")
  plot_df = data.frame(dif_df)
  #plot_df$ct1 = factor(plot_df$ct1, levels = levels(plot_df$ct1)[hcl$order])
  #plot_df$ct2 = factor(plot_df$ct2, levels = levels(plot_df$ct2)[hcl$order])
  plot_df$ct1 = factor(plot_df$ct1, levels = levels(net_names_all$Var1))
  plot_df$ct2 = factor(plot_df$ct2, levels = levels(net_names_all$Var1))
  #plot_df$ct1 = factor(plot_df$ct1, levels = rev(c("Stellate cells", "LSEC", "Endothelial cells\n(non-LSEC)",
  #                                                 "Kupffer cells", "cDCs", "Macrophages", "Cholangiocytes",
  #                                                 "pDCs", "ab-T cells", "Hepatocytes", "gd-T cells", 
  #                                                 "B cells", "Plasmablasts")))
  #plot_df$ct2 = factor(plot_df$ct2, levels = rev(c("Stellate cells", "LSEC", "Endothelial cells\n(non-LSEC)",
  #                                                 "Kupffer cells", "cDCs", "Macrophages", "Cholangiocytes",
  #                                                 "pDCs", "ab-T cells", "Hepatocytes", "gd-T cells", 
  #                                                 "B cells", "Plasmablasts")))
  pltcnts = ggplot(plot_df[order(plot_df$Freq, decreasing = F),], 
                   aes(x = ct1, y = ct2, size = abs(Freq), colour = Freq))+
    geom_point()+
    labs(x = "cell type", y = "cell type", title = paste0(n, " - healthy"))+
    guides(colour = guide_legend(title = "change in\n# interactions", 
                                 reverse = T, override.aes = list(size = 4)), 
           size = guide_none())+
    scale_colour_gradientn(colors = c("#d7191c", "#fdae61", "#ededed", "#c2a5cf", "#7b3294"), 
                           limits = c(-110, 110))+
    scale_size_continuous(range = c(0.2, 5))+
    theme_classic()+
    theme(aspect.ratio = 1,
          panel.background = element_rect(colour = "black"),
          axis.line = element_blank(),
          axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.4),
          axis.text = element_text(colour = "black", size = 8))
  print(pltcnts)
  diff_cnt_l[[n]] = pltcnts
}

diff_cnt_l$legend = cowplot::get_legend(diff_cnt_l$embolised)
diff_cnt_l$embolised = diff_cnt_l$embolised+theme(legend.position = "none")
diff_cnt_l$regenerating = diff_cnt_l$regenerating+theme(legend.position = "none")

pdf("figure_panels_rev1/fig_comm/suppFig_countsDiff.pdf", height = 8, width = 15)
print(cowplot::plot_grid(plotlist = diff_cnt_l[c(2,1,3)], ncol = 3, rel_widths = c(1,1,0.3), align = "h"))

Warning: Graphs cannot be horizontally aligned unless the axis parameter is set. Placing graphs unaligned.

dev.off()

png 
  2 

Unique interactions from E and R

print(cowplot::plot_grid(plotlist = diff_cnt_l[c(2,1,3)], 
                         ncol = 3, rel_widths = c(1,1,0.3), align = "hv"))

Warning: Graphs cannot be vertically aligned unless the axis parameter is set. Placing graphs unaligned.
Warning: Graphs cannot be horizontally aligned unless the axis parameter is set. Placing graphs unaligned.

Functions to plot interactions

# number of new unique interactions
## i.e. this ligand-receptor pair was absent at least in healthy and always includes this ct
# get interactions that are not in healthy
inter_h = inter_df$healthy$id_cp_interaction
inter_reg_sub = inter_df$regenerating[!inter_df$regenerating$id_cp_interaction %in% inter_h,
                                      1:3]
inter_emb_sub = inter_df$embolised[!inter_df$embolised$id_cp_interaction %in% inter_h, 1:3]

pve_inter = rbind(inter_reg_sub, inter_emb_sub)
pve_inter$cond = c(rep("regenerating", nrow(inter_reg_sub)),
                   rep("embolised", nrow(inter_emb_sub)))
pve_inter$mid_ct1 = match_df[pve_inter$ct1,"ct_mid"]
pve_inter$mid_ct2 = match_df[pve_inter$ct2,"ct_mid"]

# list unique interactions (all ct)
ct_inter_l = list()
for(cond in unique(pve_inter$cond)){
  ct_inter_l[[cond]] = list()
  sub_pve_inter = pve_inter[pve_inter$cond==cond,]
  sub_pve_inter = sub_pve_inter[sub_pve_inter$ct1!=sub_pve_inter$ct2,]
  
  inter_u = unique(sub_pve_inter$id_cp_interaction)
  for(i in inter_u){
    idf = sub_pve_inter[sub_pve_inter$id_cp_interaction==i, c("ct1", "ct2")]
    tab_int_nu = table(idf$ct1, idf$ct2)
    
    if(nrow(tab_int_nu)==1){
      ct_inter_l[[cond]][[i]] = rownames(tab_int_nu)
    }
    if(ncol(tab_int_nu)==1){
      ct_inter_l[[cond]][[i]] = colnames(tab_int_nu)
    }
  }
  
}
ct_inter_l = reshape2::melt(ct_inter_l)
table(ct_inter_l$value, ct_inter_l$L1)

                            
                             embolised regenerating
  ab-T cells (stress)                3            0
  activated DCs                      0            1
  B cells                            1            1
  Dividing cDCs                      1            1
  Dividing endothelial cells         3            0
  EC non-LSEC                        2            0
  IgA+ Plasma cells                  0            1
  ILC3                               0            2
  Infiltrating NK cells              2            0
  LSEC (fenestr.)                    1            1
  LSEC (interferon)                  0            1
  LSEC (remodelling)                 0            3
  Lymphatic EC                       2            0
  Monocytes (TREM2+ CD9+)            1            1
  pDCs                               4            4
  Stellate cells                    10           17
  Treg                               1            0

# list unique interactions (mid ct)
ctmid_inter_l = list()
ctmid_inter_l_ct = list()
for(cond in unique(pve_inter$cond)){
  ctmid_inter_l[[cond]] = list()
  ctmid_inter_l_ct[[cond]] = list()
  sub_pve_inter = pve_inter[pve_inter$cond==cond,]
  sub_pve_inter = sub_pve_inter[sub_pve_inter$mid_ct1!=sub_pve_inter$mid_ct2,]
  
  inter_u = unique(sub_pve_inter$id_cp_interaction)
  for(i in inter_u){
    idf = sub_pve_inter[sub_pve_inter$id_cp_interaction==i, c("mid_ct1", "mid_ct2")]
    tab_int_nu = table(idf$mid_ct1, idf$mid_ct2)
    mid_ct1_ct = sub_pve_inter[sub_pve_inter$id_cp_interaction==i, "ct1"]
    mid_ct2_ct = sub_pve_inter[sub_pve_inter$id_cp_interaction==i, "ct2"]
    
    if(nrow(tab_int_nu)==1){
      ctmid_inter_l[[cond]][[i]] = rownames(tab_int_nu)
      ctmid_inter_l_ct[[cond]][[i]] = mid_ct1_ct
    }
    if(ncol(tab_int_nu)==1){
      ctmid_inter_l[[cond]][[i]] = colnames(tab_int_nu)
      ctmid_inter_l_ct[[cond]][[i]] = mid_ct2_ct
    }
  }
  
}
ctmid_inter_l = reshape2::melt(ctmid_inter_l)
ctmid_inter_l_ct = unique(reshape2::melt(ctmid_inter_l_ct))
ctmid_inter_l = merge(ctmid_inter_l, ctmid_inter_l_ct, by = c("L2", "L1"), 
                      all = T)[,c(3,1,2,4)]
table(ctmid_inter_l$value.x, ctmid_inter_l$L1)

                
                 embolised regenerating
  B cells                1            2
  Hepatocytes            9            6
  ILC                    2            2
  LSEC                   4            3
  Mesenchymal           10           17
  other ECs              8           10
  other Mono-Mac         2            3
  pDCs                   4            4
  T cells                3            0

Plot all

filtFuncU = function(x, intdf, mutthr = 0, expthr = 0.12, minct = 3){
  coldf = x[,1:8]
  #int_ct_p_t = paste0(intdf$value, intdf$L2)
  coldf$isMain1 = apply(coldf, 1, function(y) paste0(y[2], y[1]) %in% intdf)
  coldf$isMain2 = apply(coldf, 1, function(y) paste0(y[3], y[1]) %in% intdf)
  
  # swap if isMain2
  for(i in 1:nrow(coldf)){
    if(coldf[i,"isMain2"]){
      tmpct = coldf[i,"ct2"]
      tmplr = coldf[i,"lr2"]
      tmpma = coldf[i,"isMain2"]
      
      coldf[i,"ct2"] = coldf[i,"ct1"]
      coldf[i,"lr2"] = coldf[i,"lr1"]
      coldf[i,"isMain2"] = coldf[i,"isMain1"]
      
      coldf[i,"ct1"] = tmpct
      coldf[i,"lr1"] = tmplr
      coldf[i,"isMain1"] = tmpma
    }
  }

  # interaction names
  coldf$intlr = paste0(coldf$lr1, " -\n ", coldf$lr2)
  coldf = unique(coldf[,c(2,3,11,6:10)])
  # interaction frequency
  nint = table(coldf$intlr)
  # expression and interaction filtering
  coldf = coldf[apply(coldf[,4:6], 1, function(x) any(x>=expthr)) & 
                  (coldf$intlr %in% names(nint)[nint>1] | 
                     apply(coldf[,4:6], 1, function(x) any(x>=expthr*8))),]
  # reshaping for plot
  coldf = data.table::rbindlist(list(reshape2::melt(coldf[,c(3,1,4:7)]), 
                                     reshape2::melt(coldf[,c(3,2,4:6, 8)])), use.names = F)
  coldf$variable = unlist(lapply(strsplit(as.character(coldf$variable), "_"), 
                                 function(x) x[1]))
  coldf$variable = factor(coldf$variable, levels = c("healthy","regenerating", "embolised"))
  coldf = coldf[order(coldf$value, decreasing = T),]
  coldf = coldf[!duplicated(coldf[,1:4]),]
  coldf$ct1 = gsub("Endothelial cells", "EC", coldf$ct1)
  coldf$intlr = gsub("receptor", "rec", coldf$intlr)
  coldf$intlr = gsub("complex", "comp", coldf$intlr)
  nct = table(coldf$ct1)
  coldf = coldf[coldf$ct1 %in% names(nct)[nct>=minct*3],]
  nint = table(coldf$intlr)
  coldf = coldf[coldf$intlr %in% names(nint)[nint>3],]
  
  # confirm that all interactions have at least one TRUE
  ntrue = tapply(coldf$isMain1, coldf$intlr, sum)
  coldf = coldf[coldf$intlr %in% names(ntrue)[ntrue>0]]
  
  # formatting
  coldf = coldf[order(coldf$isMain1, coldf$intlr, -coldf$value, decreasing = T),]
  coldf$ct1 = factor(coldf$ct1, levels = rev(unique(coldf$ct1)))
  coldf$intlr = factor(coldf$intlr, levels = rev(unique(coldf$intlr)))
  coldf = coldf[order(coldf$isMain1, decreasing = F),]
  
  return(coldf)
}

pltFuncU = function(x){
  plt = ggplot(x, aes(y = ct1, x = variable, fill = value, colour = isMain1))+
    facet_grid(intlr~., scales = "free_y", space = "free_y")+
    geom_tile(size = 0.4)+
    scale_x_discrete(expand = c(0,0))+
    scale_fill_gradientn(colors = gradexpcol)+
    scale_colour_manual(values = c("white", "black"))+
    guides(fill = guide_colourbar(barheight = unit(0.4, "cm")),
           colour = guide_none())+
    labs(x = "Condition", y = "Cell type", fill = "mean exp")+
    theme_classic()+
    theme(axis.text.y = element_text(colour = "black", size = 7, 
                                     angle = 0, hjust = 1, vjust = 0.5),
          axis.text.x = element_text(colour = "black", size = 6.5, 
                                     angle = 35, hjust=1, vjust=1),
          axis.title = element_text(size = 7),
          axis.line = element_blank(),
          legend.position = "bottom",
          legend.title = element_text(size = 7),
          legend.box.margin = margin(0,0,0,0),
          legend.text = element_text(size = 6),
          strip.text.y = element_text(angle = 0, size = 7, margin = margin(0,0,3,0), 
                                      hjust = 0.5, vjust = 0.5, face = "bold"),
          strip.background = element_blank())
  return(plt)
}

Plot one by one

plot_df = filtFuncU(her_allint_f[her_allint_f$id_cp_interaction %in% ct_inter_l$L2[!grepl("stress", ct_inter_l$value)],],
                    intdf = paste0(ct_inter_l$value, ct_inter_l$L2), 
                    expthr = 0.2, minct = 1)

Using intlr, ct1, isMain1 as id variables
Using intlr, ct2, isMain2 as id variables

imm_plt = pltFuncU(plot_df)
pdf("figure_panels_rev1/fig_comm/mainFig_InteractUnique.pdf", 
    height = 20, width = 2.8)
print(imm_plt)
dev.off()

null device 
          1 

plot_df = filtFuncU(her_allint_f[her_allint_f$id_cp_interaction %in% ctmid_inter_l$L2[!grepl("stress", ctmid_inter_l$value.y)],],
                   intdf = paste0(ctmid_inter_l$value.y, ctmid_inter_l$L2), 
                   expthr = 0.2, minct = 3)

Using intlr, ct1, isMain1 as id variables
Using intlr, ct2, isMain2 as id variables

imm_plt = pltFuncU(plot_df)
pdf("figure_panels_rev1/fig_comm/mainFig_InteractUnique_mid.pdf", 
    height = 32, width = 2.8)

print(imm_plt)
dev.off()

null device 
          1 

Plot one by one, all in H vs PVE network

plot_df = filtFuncU(her_allint_f[her_allint_f$id_cp_interaction %in% ct_inter_l$L2[!grepl("stress", ct_inter_l$value)],],
                    intdf = paste0(ct_inter_l$value, ct_inter_l$L2), 
                    expthr = 0.2, minct = 3)

Using intlr, ct1, isMain1 as id variables
Using intlr, ct2, isMain2 as id variables

folder = "figure_panels_rev1/fig_comm/heatmaps_interact_indiv/"
for(i in 1:length(unique(plot_df$intlr))){
  iii = unique(plot_df$intlr)[i]
  sub_plot_df = plot_df[plot_df$intlr==iii,]
  sub_plot_df = sub_plot_df[order(sub_plot_df$isMain1, sub_plot_df$value, decreasing = T),]
  sub_plot_df$ct1 = factor(sub_plot_df$ct1, levels = rev(unique(sub_plot_df$ct1)))
  sub_plot_df = sub_plot_df[order(sub_plot_df$isMain1, decreasing = F),]
  
  pdf(paste0(folder, "mainFig_InteractUnique_", i, ".pdf"), height = 3, width = 4)
  print(pltFuncU(sub_plot_df))
  dev.off()
}

plot_df = filtFuncU(her_allint_f[her_allint_f$id_cp_interaction %in% ctmid_inter_l$L2[!grepl("stress", ctmid_inter_l$value.y)],],
                   intdf = paste0(ctmid_inter_l$value.y, ctmid_inter_l$L2), 
                   expthr = 0.2, minct = 3)

Using intlr, ct1, isMain1 as id variables
Using intlr, ct2, isMain2 as id variables

folder = "figure_panels_rev1/fig_comm/heatmaps_interact_indiv/"
for(i in 1:length(unique(plot_df$intlr))){
  iii = unique(plot_df$intlr)[i]
  sub_plot_df = plot_df[plot_df$intlr==iii,]
  sub_plot_df = sub_plot_df[order(sub_plot_df$isMain1, sub_plot_df$value, decreasing = T),]
  sub_plot_df$ct1 = factor(sub_plot_df$ct1, levels = rev(unique(sub_plot_df$ct1)))
  sub_plot_df = sub_plot_df[order(sub_plot_df$isMain1, decreasing = F),]
  
  pdf(paste0(folder, "mainFig_InteractUnique_mid_", i, ".pdf"), height = 3, width = 2.8)
  print(pltFuncU(sub_plot_df))
  dev.off()
}

Plot changes in number of interactions per condition - mid resolution

plot_df = filtFunc(her_allint_f[her_allint_f$id_cp_interaction %in% unique(edge_min$id_cp_interaction),],
                   expthr = 0.15, minct = 2)

Using intlr, ct1 as id variables
Using intlr, ct2 as id variables

Panels cirrhosis

UMAP Immune cells

UMAP Mono/Mac cells

m_df = readRDS(file = "results/cirrhosis/umap_m_df.RDS")

cols = MetBrewer::met.brewer("Klimt", length(unique(m_df$all_annot)))
names(cols) = unique(m_df$all_annot)

m_df = m_df[sample(1:nrow(m_df), nrow(m_df), replace = F),]
plt = ggplot()+
  geom_point(aes(x = UMAP_1, y = UMAP_2, colour = all_annot), 
             m_df, size = 0.2)+
  guides(colour = guide_legend(override.aes = list(size = 3)))+
  scale_colour_manual(values = cols[order(names(cols))])+
  theme_void()+
  theme(aspect.ratio = 1,
        legend.text = element_text(size = 6),
        legend.key.size = unit(0.4, "cm"),
        legend.title = element_blank())

pdf("figure_panels_rev1/fig_cir/all_mono_umap.pdf", height = 2.2, width = 4.2)
print(plt)
dev.off()

UMAP T/NK cells

t_df = readRDS(file = "results/cirrhosis/umap_t_df.RDS")

cols = MetBrewer::met.brewer("Juarez", length(unique(t_df$all_annot)))
names(cols) = unique(t_df$all_annot)

t_df = t_df[sample(1:nrow(t_df), nrow(t_df), replace = F),]
plt = ggplot()+
  geom_point(aes(x = UMAP_1, y = UMAP_2, colour = all_annot), 
             t_df, size = 0.2)+
  guides(colour = guide_legend(override.aes = list(size = 3)))+
  scale_colour_manual(values = cols[order(names(cols))])+
  theme_void()+
  theme(aspect.ratio = 1,
        legend.text = element_text(size = 6),
        legend.key.size = unit(0.375, "cm"),
        legend.title = element_blank())

pdf("figure_panels_rev1/fig_cir/all_tnk_umap.pdf", height = 2.2, width = 4.2)
print(plt)
dev.off()

Markers myeloid cells

mat_m = readRDS(file = "results/cirrhosis/mat_myeloid_markers.RDS")

cols_pal = colorRampPalette(rev(brewer.pal(n = 9, name = "RdBu")))(100)

pdf("figure_panels_rev1/fig_cir/heat_mk_mye.pdf", height = 2.2, width = 6.3)
pheatmap::pheatmap(mat_m, clustering_method = "ward.D", color = cols_pal,
                   treeheight_col = 0, treeheight_row = 20, fontsize = 6.3)
dev.off()

Markers lymphoid cells

mat_l = readRDS("results/cirrhosis/mat_lymphoid_markers.RDS")

cols_pal = colorRampPalette(rev(brewer.pal(n = 9, name = "RdBu")))(100)

pdf("figure_panels_rev1/fig_cir/heat_mk_lym.pdf", height = 2.2, width = 6.3)
pheatmap::pheatmap(mat_l, clustering_method = "ward.D", color = cols_pal,
                   treeheight_col = 0, treeheight_row = 20, fontsize = 6.3)
dev.off()

Proportions immune cells

imm_props = readRDS(file = "results/cirrhosis/imm_props_dat.RDS")
imm_pvals = readRDS("results/cirrhosis/imm_props_pval.RDS")

plt = ggplot(imm_props, aes(x = Condition, y = Freq, group = Condition, colour = Condition))+
  facet_wrap(~Var1, scales = "free")+
  geom_jitter(position = position_jitterdodge(jitter.width = 0.3, dodge.width = 1))+
  stat_summary(fun.data = mean_se, position = position_dodge(width = 1), 
               alpha = 0.35, colour = "black")+
  theme_bw()+
  theme(legend.position = "none")


colnames(imm_pvals)[1] = "embolized"

imm_props$Condition[imm_props$Condition=="embolised"] = "embolized"
imm_props$Condition = factor(imm_props$Condition, 
                             levels = c("healthy", "regenerating", "embolized"))

pval_pos = imm_props %>%
  group_by(Condition, Var1) %>%
  summarise_at(vars("Freq"), mean)
pval_pos$posy = pval_pos$Freq*1.02
pval_pos$posx = as.numeric(pval_pos$Condition)+0.3
pval_pos = pval_pos[pval_pos$Condition!="healthy",]
pval_pos$pval = diag(as.matrix(imm_pvals[pval_pos$Var1,as.character(pval_pos$Condition)]))
pval_pos$issig = ifelse(pval_pos$pval<=0.05, "*", "")

plt = ggplot()+
  facet_wrap(~Var1, scales = "free")+
  geom_jitter(data = imm_props, mapping = aes(x = Condition, y = Freq, 
                                              group = Condition, colour = Condition),
              position = position_jitterdodge(jitter.width = 0.3, dodge.width = 1), size = 0.5)+
  stat_summary(data = imm_props, mapping = aes(x = Condition, y = Freq, 
                                              group = Condition, colour = Condition),
               fun.data = mean_se, position = position_dodge(width = 1), 
               alpha = 0.6, colour = "black", size = 0.25)+
  geom_text(data = pval_pos, mapping = aes(x = posx, y = posy, label = issig), size = 4)+
  scale_colour_manual(values = colcond, drop = T, limits = names(colcond)[c(1,2,4)])+
  scale_x_discrete(labels = c("H", "R", "E"))+
  labs(y = "proportion % (immune)")+
  theme_bw()+
  theme(legend.position = "none",
        axis.title = element_text(size = 6, colour = "black"),
        axis.text.x = element_text(size = 6, colour = "black"),
        axis.text.y = element_text(size = 6, colour = "black"),
        strip.text = element_text(size = 6),
        strip.background = element_rect(fill = "white"),
        panel.grid.minor = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.major.y = element_line(size = 0.3, colour = "grey70"))

#pdf("figure_panels_rev1/fig_cir/immune_all_prop.pdf", height = 7.5, width = 11.5)
print(plt)
#dev.off()

pdf("figure_panels_rev1/fig_cir/immune_all_prop_resized.pdf", height = 5.5, width = 7)
print(plt)
dev.off()

Proportions endothelial cells

plt = ggplot()+
  facet_wrap(~Var1, scales = "free")+
  geom_jitter(data = imm_props, mapping = aes(x = Condition, y = Freq, 
                                              group = Condition, colour = Condition),
              position = position_jitterdodge(jitter.width = 0.3, dodge.width = 1), size = 0.5)+
  stat_summary(data = imm_props, mapping = aes(x = Condition, y = Freq, 
                                              group = Condition, colour = Condition),
               fun.data = mean_se, position = position_dodge(width = 1), 
               alpha = 0.6, colour = "black", size = 0.25)+
  geom_text(data = pval_pos, mapping = aes(x = posx, y = posy, label = issig), size = 4)+
  scale_colour_manual(values = colcond, drop = T, limits = names(colcond)[c(1,2,4)])+
  scale_x_discrete(labels = c("H", "R", "E"))+
  labs(y = "proportion % (immune)")+
  theme_bw()+
  theme(legend.position = "none",
        axis.title = element_text(size = 6, colour = "black"),
        axis.text.x = element_text(size = 6, colour = "black"),
        axis.text.y = element_text(size = 6, colour = "black"),
        strip.text = element_text(size = 6),
        strip.background = element_rect(fill = "white"),
        panel.grid.minor = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.major.y = element_line(size = 0.3, colour = "grey70"))

#pdf("figure_panels_rev1/fig_cir/immune_all_prop.pdf", height = 7.5, width = 11.5)
print(plt)
#dev.off()

pdf("figure_panels_rev1/fig_cir/immune_all_prop_resized.pdf", height = 5.5, width = 7)
print(plt)
dev.off()

png 
  2 

MP signatures

mp_sig_df = readRDS(file = "results/cirrhosis/mono_cirr_ct_signatures.RDS")

plot_df = reshape2::melt(mp_sig_df)
plot_df$Condition[plot_df$Condition=="embolised"] = "embolized"
plot_df$Condition = factor(plot_df$Condition, 
                           levels = c("healthy","regenerating","embolized"))
sub_plot_df = plot_df[plot_df$variable %in% c("cirr_ct_MPs_4","cirr_ct_MPs_5"),]
sub_plot_df$variable = gsub("cirr_ct_MPs_4", "SAMac (1)", sub_plot_df$variable)
sub_plot_df$variable = gsub("cirr_ct_MPs_5", "SAMac (2)", sub_plot_df$variable)

plt = ggplot(sub_plot_df, aes(x = mono_annot, y = value, fill = Condition))+
  facet_grid(variable~mono_annot, scales = "free")+
  geom_hline(yintercept = 0, linetype = "dashed", size = 0.3)+
  geom_violin(scale = "count")+
  scale_fill_manual(values = colcond, drop = T, limits = names(colcond)[c(1,2,4)])+
  theme_classic()+
  theme(axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.text.y = element_text(size = 6, colour = "black"),
        axis.title.y = element_text(size = 6.5),
        strip.text = element_text(size = 7),
        panel.border = element_rect(colour = "black", fill = NA),
        axis.line = element_blank(),
        legend.position = c(0.048, 0.865),
        legend.background = element_blank(),
        legend.key.size = unit(0.3, "cm"),
        legend.text = element_text(size = 6.5),
        legend.title = element_text(size = 7))

pdf("figure_panels_rev1/fig_cir/mp_sig_vio.pdf", height = 3, width = 10.35)
print(plt)
dev.off()

Endothelial signatures

mp_sig_df = readRDS(file = "results/cirrhosis/mono_cirr_ct_signatures.RDS")

plot_df = reshape2::melt(mp_sig_df)

Using Condition, mono_annot as id variables

plot_df$Condition[plot_df$Condition=="embolised"] = "embolized"
plot_df$Condition = factor(plot_df$Condition, 
                           levels = c("healthy","regenerating","embolized"))
sub_plot_df = plot_df[plot_df$variable %in% c("cirr_ct_MPs_4","cirr_ct_MPs_5"),]
sub_plot_df$variable = gsub("cirr_ct_MPs_4", "SAMac (1)", sub_plot_df$variable)
sub_plot_df$variable = gsub("cirr_ct_MPs_5", "SAMac (2)", sub_plot_df$variable)

plt = ggplot(sub_plot_df, aes(x = mono_annot, y = value, fill = Condition))+
  facet_grid(variable~mono_annot, scales = "free")+
  geom_hline(yintercept = 0, linetype = "dashed", size = 0.3)+
  geom_violin(scale = "count")+
  scale_fill_manual(values = colcond, drop = T, limits = names(colcond)[c(1,2,4)])+
  theme_classic()+
  theme(axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.text.y = element_text(size = 6, colour = "black"),
        axis.title.y = element_text(size = 6.5),
        strip.text = element_text(size = 7),
        panel.border = element_rect(colour = "black", fill = NA),
        axis.line = element_blank(),
        legend.position = c(0.048, 0.865),
        legend.background = element_blank(),
        legend.key.size = unit(0.3, "cm"),
        legend.text = element_text(size = 6.5),
        legend.title = element_text(size = 7))

pdf("figure_panels_rev1/fig_cir/mp_sig_vio.pdf", height = 3, width = 10.35)
print(plt)
dev.off()

null device 
          1 

Hepatocyte signatures

end_sig_df = readRDS(file = "results/cirrhosis/endo_cirr_ct_signatures.RDS")

plot_df = reshape2::melt(end_sig_df)

Using Condition, endo_simp as id variables

plot_df$Condition[plot_df$Condition=="embolised"] = "embolized"
plot_df$Condition[plot_df$Condition=="embolised"] = "embolized"
plot_df$Condition = factor(plot_df$Condition, 
                           levels = c("healthy","regenerating","embolized"))
sub_plot_df = plot_df[plot_df$variable %in% c("cirr_ct_Endothelia_6","cirr_ct_Endothelia_7"),]
sub_plot_df$variable = gsub("cirr_ct_Endothelia_6", 
                            "SAEndothelial (1)", sub_plot_df$variable)
sub_plot_df$variable = gsub("cirr_ct_Endothelia_7", "SAEndothelial (2)",
                            sub_plot_df$variable)

plt = ggplot(sub_plot_df, aes(x = endo_simp, y = value, fill = Condition))+
  facet_grid(variable~endo_simp, scales = "free")+
  geom_hline(yintercept = 0, linetype = "dashed", size = 0.3)+
  geom_violin(scale = "count")+
  scale_fill_manual(values = colcond, drop = T, limits = names(colcond)[c(1,2,4)])+
  theme_classic()+
  theme(axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.text.y = element_text(size = 6, colour = "black"),
        axis.title.y = element_text(size = 6.5),
        strip.text = element_text(size = 7),
        panel.border = element_rect(colour = "black", fill = NA),
        axis.line = element_blank(),
        legend.position = "bottom",
        legend.background = element_blank(),
        legend.margin = margin(0,0,0,0),
        legend.key.size = unit(0.3, "cm"),
        legend.text = element_text(size = 6.5),
        legend.title = element_text(size = 7))

pdf("figure_panels_rev1/fig_cir/endo_sig_vio.pdf", height = 2.8, width = 11)
print(plt)
dev.off()

null device 
          1